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" PROBIT":  A  PROBIT  ANALYSIS  PROGRAM  FOR  THE  DRES  COMPUTER  FACILITY 


by 

J.  McFee,  C.E.  Mendoza, 
J.  Smith  and  D.  Rigby 


ABSTRACT 


This  report  describes  an  interactive  computer  program,  PROBIT,  which 
performs  automated  probit  analysis.  The  code  is  written  in  the  FORTRAN  IV 
language  using  a  modified  iterative  maximum  likelihood  algorithm  and  runs  on  the 
DRES  Honeywell  DPS-8/70  and  VAX  11/780  computers. 

A  series  of  tests  were  performed  to  compare  PROBIT  with  other  probit 
analysis  programs.  These  include  the  old  DRES  PROBT  (modified  I8M  PROBT)  and 
programs  contained  in  the  commerical  statistical  packages  S103  and  SAS.  PROBIT 
clearly  outperformed  DRES  PROBT  and  performed  as  well  as  S103  and  as  well  or 
slightly  better  than  SAS. 

The  input  methods  and  output  format  of  PROBIT  can  be  modified  whereas  those 
of  most  commerical  statistical  packages  cannot.  Finally,  PR08IT  cost 
substantially  less  to  develop  than  the  purchase  price  of  typical  commerical 
statistical  packages  containing  a  good  probit  analysis  program. 
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"PROBIT":  A  PROBIT  ANALYSIS  PROGRAM  FOR  THE  ORES  COMPUTER  FACILITY 


by 

J.  McFee,  C.E.  Mendoza, 
J.  Smith  and  D.  Rigby 


1.0  INTRODUCTION 


There  are  numerous  ways  to  quantify  the  effect  of  a  substance  on  an 
organism.  One  of  the  most  widely  used  Is  the  ED50,  the  quantity  of  substance 
which  produces  the  effect  of  Interest  In  50%  of  the  population. 

Most  experimental  designs  used  to  determine  the  EDS0  consist  of  a  finite 
number  of  treatment  groups,  each  of  which  comprises  a  finite  number  of  organisms. 
Orlganlsms  within  each  group  are  given  the  same  quantity  of  substance  and  the 
quantity  differs  from  group  to  group.  Usually  the  number  of  treatment  groups  and 
the  number  of  organisms  within  a  group  are  small  U10). 

A  variety  of  methods  exists  to  estimate  ED50  levels  from  the  data.  These 
Include  the  methods  of  extreme  lethal  doses,  the  Dragstedt-Behrens  method, 
Reed-Muench,  Spearman-Karber,  moving  average,  and  probit  analysis  (Finney  1971). 
Of  these,  only  probit  analysis  Is  completely  self  consistent.  The  other  methods 
will  commonly  give  the  same  answer  as  probit  analysis,  but  lead  to  serious  errors 
If  the  dose  range  Is  very  asymmetric  about  ln(EDS0) . 
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There  are  several  methods  for  performing  probit  analysis.  Some  are 
graphical  and  therefore  not  readily  suited  for  automated  data  analysis.  Others 
cannot  estimate  ED  levels  other  than  the  EDSl,  or  do  not  determine  uncertainties 
of  estimated  quantities.  A  complete  discussion  of  the  various  ways  of  estimating 
ED  levels  may  be  found  in  Finney  (1971). 

In  the  past,  ORES  carried  out  probit  analysis  using  a  program,  PROBT 
(Chenier  et  al.  1962),  based  on  a  modified  form  of  the  subroutine  PROBT  (IBM 
PROBT)1  which  is  contained  in  the  IBM  Scientific  Subroutine  Library'-.  The 
subroutine  is  based  on  a  modified  version  of  the  iterative  maximum  likelihood 
method  of  Finney  (1971).  This  is  a  popular  method  for  automated  probit  analysis 
since  it  requires  no  further  user  intervention  once  data  has  been  provided  and  is 
quite  robust.  (Robustness  of  an  algorithm  or  program  is  a  measure  of  its  ability 
to  converge  rapidly  to  a  solution,  and  to  function  without  producing  meaningless 
numbers  or  encountering  fatal  errors.)  A  number  of  faults  had  been  demonstrated 
with  the  DRES  PROBT  program.  ORES  PROBT  did  not  function  well  with  data  sets 
which  included  a  number  of  null  and/or  100%  response  rates.  Further,  the  program 
calculated  only  the  EDbl);  whereas  in  some  applications,  other  ED  levels  such  as 
the  ED p  E995  or  EU99  were  often  required.  Recently,  it  was  also  found  that  the 
code  contains  errors  in  calculation.  These  included  incorrect  estimates  of  the 
number  of  degrees  of  freedom  and  the  variance  of  the  regression  slope,  as  well  as 
some  generally  incorrect  equations  which  fortuituously  yielded  correct  answers 
only  for  ED5Q.  An  even  more  serious  shortcoming  of  DRES  PROBT  was  the  fact  that 
it  Issued  no  warnings  for  a  number  of  cases  where  results  would  be  suspect,  such 
as  g  >  1  or  heterogeneous  data.  These  shortcomings,  together  with  the  fact  that 
the  original  program  had  been  written  for  the  obsolete  IBM  1130  computer  in  an 
outmoded  programming  language  (FORTRAN  II),  prompted  us  to  design  a  new  program 
that  would  eliminate  the  above  problems,  be  cost  effective,  run  on  the  DRES 
Honeywell  DPS-8  computer  and  be  transportable  to  different  computers. 


1  In  this  report,  the  modified  IBM  PROBT  will  be  referred  to  as  DRES  PROBT. 

2  System/360  Scientific  Subroutine  Package  (360A-CM-03X)  Version  3  Programmers 

Manual,  IBM  Technical  Publications  Department,  White  Plains,  N.Y.,  44  (1968). 
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The  method  of  probit  analysis  used  for  the  new  program  Is  also  based  on  the 
iterative  maximum  likelihood  method  of  Finney.  Robustness  for  Finney's  method  is 
determined  primarily  by  two  key  items  in  the  program.  The  first  is  the  method  of 
handling  null  and  100%  responses  in  the  initial  iteration.  Null  responses  must 
be  eliminated  for  the  first  iteration  since  they  correspond  to  probits  equal  to 
minus  infinity.  The  inclusion  of  100%  responses  can  also  cause  instabilities  in 
the  program  particularly  if  several  null  and  100%  responses  are  present  in  the 
data,  but  it  is  debatable  whether  they  should  be  eliminated  from  the  initial 
iteration  or  merely  weighted  differently  than  the  other  responses.  The  second 
key  item  is  the  method  of  testing  for  convergence.  There  are  several  quantities 
that  can  be  used  to  test  convergence  but  which  one  provides  the  most  accurate 
answers  for  the  minimum  number  of  iterations  is  not  clear.  Since  most  successful 
commerical  probit  analysis  programs  have  their  own  methods  of  handling  these  two 
operations  and  they  do  not  provide  details  of  their  algorithms,  one  can  only 
determine  how  to  handle  convergence  and  null/100%  responses  by  trying  different 
approaches.  After  extensive  tests,  it  was  decided  to  weigh  both  the  null  and 
100%  responses  by  zero  in  the  first  iteration  and  to  use  the  sensitivity  (slope 
of  the  regression  line)  as  the  test  of  convergence. 


The  new  PROBIT  code  is  an  interactive  FORTRAN  IV  program  which  yields  a 
complete  analysis  of  probit  experiments,  or,  if  desired,  a  condensed  summary  of 
results.  The  program  provides  a  far  more  thorough  analysis  and  is  markedly  less 
sensitive  to  the  dynamic  range  of  the  input  data  than  many  personal  computer 
probit  programs  (Lleberman,  1983)  and  IBM  PROBT. 

There  are  equally  sophisticated  programs  available  for  mainframe  computers, 
but  these  were  not  deemed  cost  effective  for  use  at  DRES,  since  they  are  a  part 
of  large  statistical  program  packages,  such  as  SAS3  and  S103u,  which 

3  Statistical  Analysis  System,  SAS  Institute  Inc.,  SAS  Circle,  Box  8000,  Cary, 
N.C. 

u  S103  Statistical  Package,  Engineering  and  Statistical  Research  Institute, 
Agriculture  Canada,  Ottawa,  K1A  0C6 
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must  be  purchased  or  rented  in  their  entirety.  Also,  these  programs  cannot 
normally  be  modified.  This  is  a  significant  problem,  since  it  was  found  that  a 
number  of  individuals  at  DRES  wanted  certain  output  formats  which  these  codes  did 
not  allow.  This  would  be  a  conti ning  problem  as  scientists  and  their  tastes 
change  with  time.  A  DRES  program,  however,  could  be  tailored  to  suit  needs  now 
and  be  modified  later. 

An  overview  of  the  mainline  program  is  discussed  in  Section  2.0. 
Sections  2.1  and  2.2  discuss  the  pertinent  input  and  output  parameters. 
Section  3  describes  the  subprograms  required  by  the  main  program.  Results  of 
tests  of  the  program,  including  a  comparison  of  a  number  of  probit  analysis 
programs  are  presented  in  Section  4.  The  underlying  theory  behind  the  program  is 
described  in  Appendices  1  and  2.  A  description  of  operating  instructions  for  the 
program  as  implemented  on  the  DRES  Honeywell  DPS-8/70  and  Ordnance  Detection 
Group  (ODG)  VAX  11/780  is  given  in  Appendix  3.  A  sample  input  data  file  and 
sample  output  are  given  in  Appendix  4.  Finally,  Appendix  5  contains  a  listing  of 
PROBIT. 

I 

2.0  MAIN  PROGRAM 


Figure  1  shows  the  hierarchial  chart  of  the  probit  analysis  program. 
Figure  2  is  the  flowchart  of  the  main  program,  PROBIT.  The  main  program  controls 
the  input  of  data,  the  output  of  results  and  the  actual  iterative  maximum 
likelihood  calculation.  Once  the  necessary  data  have  been  entered,  the  iterative 
procedures  begins.  It  continues  until  convergence  occurs  (Appendix  1)  or  until 
the  number  of  iterations  exceeds  a  specified  value,  indicating  nonconvergence. 

The  user  has  the  option  of  printing  intermediate  results  of  the  iterative 
regression  at  the  end  of  each  iteration.  This  option  is  controlled  by  a  switch 
parameter  during  the  input  portion  of  the  program,  and  is  particularly  useful  in 
analysing  the  results  obtained  from  data  which  produced  nonconvergence  or  results 
which  should  be  viewed  with  scepticism. 

The  output  comprises  some  of  the  best  features  of  a  number  of  probit 
analysis  programs  designed  for  mainframe  computers,  notably  SAS  and  S103. 
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2.1  Input 

All  of  the  input  required  for  the  probit  analysis  program  is  "interactive" 
and  is  entered  directly  into  the  computer,  by  means  of  the  terminal  keyboard. 
The  data  is  entered  in  "free"  format  (excluding  the  alphanumeric  strings),  with 
the  numbers  separated  by  a  comma  or  a  blank  space.  The  program  itself  need  not 
be  altered  for  execution;  and  the  user  is  told  what  values  to  enter  with  default 
values  provided  where  they  might  be  needed.  If  the  input  is  entered  incorrectly, 
execution  may  be  interrupted  and  the  values  re-entered.  Some  initial  data  are 
printed  with  the  final  results.  These  are  shown  in  Table  1.  (In  Tables  1  or  2, 
symbols  in  brackets  refer  to  the  corresponding  variables  used  in  Appendices  1  and 
2.) 


2.2  Output 


The  output  from  a  sample  analysis  is  given  in  Appendix  4.  Intermediate 
results  have  been  printed  by  setting  the  IYNSW  switch  to  YES.  The  first  three 
iterations  are  shown  for  brevity. 

The  output  for  the  program  PROBIT  consists  of  an  initial  list  of  the  input 
data,  followed  by  intermediate  results  if  required  for  further  analysis,  and 
finally  the  results  of  the  analysis  (see  Appendices  1  and  2).  These  include  the 
doses,  dose  metameters,  number  of  responses,  probits,  NEDS,  standard  errors,  and 
weights  used  for  each  trial;  the  means  and  fiducial  limits  of  the  estimated 
responses,  the  sensitivity  with  its  standard  error;  the  value  of  chi-square  (x2) 
for  the  corresponding  number  of  degrees  of  freedom;  the  heterogeneity  factor;  the 
values  of  6  and  T  (see  below);  the  effective  dose  (ED)  levels;  and  finally,  the 
upper  and  lower  fiducial  limits  corresponding  to  the  given  ED  tolerance  levels. 
All  data  are  displayed  in  table  form  with  suitable  headings,  and  the  title  is 
printed  at  the  top  of  each  new  page.  Error  messages  have  also  been  included 
where  they  may  be  needed.  The  final  output  quantities  are  summarized  in 
Table  2. 
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Data  which  may  be  used  to  calculate  relative  potency  are  printed  in  a  format 
suitable  for  direct  entry  to  a  program,  RPOTC  (Appendix  6).  The  theory  behind 
the  program  is  beyond  the  scope  of  this  report,  but  is  found  in  Chenier,  1966  and 
Finney,  1971. 


3.0  FUNCTION  SUBPROGRAMS 

Four  function  subprograms  (one  of  which  calls  two  more  function  subprograms) 
are  required  by  the  main  program  to  transform  between  NED  values  and  probability 
of  response,  and  to  compute  values  of  the  x2  and  Student-t  distributions.  The 
subprograms  are  briefly  described  below. 

3.1  Name:  Function  XNED(Q)  (Figure  3) 

Purpose:  The  function  XNED(Q)  is  called  by  the  main  program  to  calculate 
the  normal  equivalent  deviate  (NED)  value  for  a  probability  of  no  response,  Q. 
XNED(Q)  has  a  maximum  error  of  ±0.0004  over  the  given  range  (0.0  <  Q  <  1.0). 
When  Q  is  outside  of  this  range,  values  are  set  accordingly:  i.e.,  if  Q  =  0, 
NED  =  5.0;  if  Q  =  1,  NED  =  -5.0.  The  transformation  used  for  this  function  is 
from  Hastings  (1955). 

Input  Parameters:  The  one  input  parameter  for  this  function  is  passed 
directly  from  the  mainline  program. 

0  -  real;  probability  of  no  response 

Output  Parameters:  None  (except  possible  error  messages). 

3.2  Name:  Function  RESP(Y)  (Figure  4) 

Purpose:  This  function  is  called  from  the  mainline  program  to  calculate  the 
probability  of  response,  P,  given  the  NED  value,  Y.  RESP(Y)  has  a  maximum  error 
of  ±3  x  10~7  and  uses  the  transformation  from  Hastings  (1955). 
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Input  Parameters:  The  only  input  parameter  for  the  function  RESP(Y)  is 
passed  from  the  main  program. 


Y  -  real;  normal  equivalent  deviate  value  (NED) 


Output  Parameters:  None . 


3.3  Name:  Function  PCHISQ(CHISQR,  NFREE)  (Figure  5) 


Purpose:  This  function  is  called  from  the  main  program,  to  evaluate  the 
probability  of  exceeding  x2,  given  the  estimate  of  reduced  x2»  CHISQR,  and  the 
number  of  degrees  of  freedom,  NFREE  (Bevington,  1969).  It  should  be  noted  that 
for  NFREE  odd  and  x2  >50,  the  calculation  is  only  an  approximation. 


Input  Parameters:  The  following  parameters  are  passed  from  the  mainline 


program. 


CHISQR  -  real;  estimate  of  reduced  x2 
NFREE  -  integer;  number  of  degrees  of  freedom 


Output  Parameters:  None,  except  possible  error  messages. 


3.4  Name:  Function  GAMMS(X)  (Figure  6) 


Purpose:  The  function  GAMMS(X)  is  called  by  PCHISQ  to  evaluate  the  gamma 
function  for  integers  and  half-integers  (Bevington,  1969). 


Input  Parameters:  The  only  parameter  for  this  function  is  passed  from  the 
function  PCHISQ. 


X  -  real;  real  form  of  the  integer  or  half-integer  argument. 


Output  Parameters:  None. 
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3.5  Name:  Function  FACTOR(N)  (Figure  7) 


Purpose:  This  function  subprogram  is  called  by  GAMMS  to  calculate  the 
factorial  function  for  integers  (Bevington,  1969).  This  is  used  when  the 
argument  is  Integral,  and  thus  when  the  gamma  function  is  identical  to  the 
factorial  function:  FACTOR(N)  =  N!.  For  N  less  than  11,  this  computation  is 
straightforward  multiplication,  but  logarithms  are  needed  for  larger  values  in 
order  to  preserve  precision. 


Input  Parameters:  The  only  input  parameter  used  in  this  function  is  passed 
from  GAMMS ( X ) . 


N  -  integer;  integer  argument  used  for  calculations. 


Output  Parameters:  None. 


3.6  Name:  Function  TPROB (FREE .ALPHA)  (Figure  8) 


Purpose:  This  function  is  called  from  the  main  program  to  find  the  value  of 
1 T '  from  a  Student' s-t  distribution  of  FREE  degrees  of  freedom,  vrfiich  is  exceeded 
with  probability  ALPHA  (2  sided  test).  The  T  distribution  table  (Thompson,  1941) 
is  stored  in  a  data  file,  and  is  read  as  input  into  a  2-dimensional  array.  The 
values  are  then  read  and  interpolation  is  carried  out  where  necessary  to  get  an 
accurate  value  of  'T * .  Tests  are  conducted  throughout  the  program  to  determine 
special  cases  for  FREE  and  ALPHA.  If  FREE  «;  1  or  ALPHA  <;  0.5,  the  program 

terminates.  If  FREE  >  120  the  value  of  FREE  =  ®  is  used,  and  if  ALPHA  ,*  0.995 

the  value  of  ALPHA  =  0.995  is  used.  Tests  are  also  used  to  determine  if 

interpolation  is  necessary,  and  if  so,  whether  it  is  needed  for  either  the  ALPHA 
or  FREE  variables,  or  both.  A  simple  first  order  (linear)  interpolation  is 


used. 


Input  Parameters:  As  well  as  the  two  parameters  passed  from  the  main 
program  there  is  the  data  array  read  from  a  data  file.  Thus,  the  following  are 
all  classed  as  input: 
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FREE  -  real ;  number  of  degrees  of  freedom 
ALPHA  -  real;  confidence  level  of  fiducial  limits 
T  -  real;  data  array  of  *T *  distribution  table 

Output  Parameters:  None,  except  possible  error  messages. 


4.0  TESTS  OF  THE  PR0B1T  PROGRAM 


In  order  to  test  how  well  the  new  probit  analysis  program  performed,  sample 
data  were  analysed  using  the  PROBIT,  DRES  PROBT,  SAS  and  S103  programs.  The 
first  two  programs  are  based  on  modified  forms  of  the  iterative  maximum 
likelihood  method  of  Finney  (1971).  No  information  is  available  on  the  specific 
algorithms  used  for  SAS  and  S103,  since  most  commerical  statistical  packages 
heavily  guard  their  methods.  However,  because  of  the  popularity  of  Finney's 
method,  it  is  suspected  that  they  too  employ  a  form  of  it. 

4.1  Comparison  of  "PROBIT"  and  "DRES  PROBT" 

The  first  test  compared  the  DRES  programs,  PROBIT  and  DRES  PROBT. 

The  data  sets  which  were  used  are  listed  in  Table  3.  The  data  sets  are 
typical  of  those  analysed  at  DRES  and  provide  a  mixture  of  homogeneous  and 
heterogeneous  data. 

Estimates  of  representative  quantities  estimated  by  both  programs  are 
presented  in  Table  4.  Other  quantities,  such  as  other  ED  levels  and 
corresponding  fiducial  limits,  are  available  from  the  new  program  PROBIT  but 
those  listed  are  sufficient  to  allow  a  satisfactory  comparison  of  the  programs. 

Both  programs  consistently  agree  on  only  two  quantities,  sensitivity  and 
EDS(,  estimates.  It  is  immediately  clear  that  DRES  PROBT  miscalculates  the  number 
of  degrees  of  freedom,  j,  since  this  should  be  the  number  of  samples  in  a  series 
minus  2.  The  reason  for  the  miscalculation  is  understood,  but  a  discussion  of  it 
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is  beyond  the  scope  of  this  report.  Also,  DRES  PROBT  does  not  estimate  "g"  when 
data  are  found  to  be  homogeneous,  nor  does  it  estimate  other  ED  values  or  their 
fiducial  limits. 


The  two  programs  disagreed  on  their  estimates  of  reduced  x2*  This  would  be 
expected  for  the  series  where  DRES  PROBT  miscalculated  v,  but  surprisingly 
occurred  even  when  v  was  correct.  The  estimates  of  the  standard  error  in  the 
sensitivities  agreed  in  all  series  except  for  the  first,  where  the  DRES  PROBT 

extimate  was  1.96  times  larger  than  that  of  PROBIT.  This  can  be  explained, 
however,  by  noting  that  in  the  first  series,  DRES  PROBT  declared  the  data  to  be 
heterogeneous  whereas  PROBIT  determined  it  to  be  homogeneous.  Thus,  DRES  PROBT 
multiplied  the  PROBIT  estimate  of  the  standard  error  by  the  square  root  of  the 
heterogeneity  factor  (=1.96).  The  discrepancy  in  heterogeneity  related  to  the 

error  in  DRES  PROBT  caused  the  number  of  degrees  of  freedom  to  be  miscalculated. 

Finally,  PROBIT  generally  converged  much  more  quickly  than  DRES  PROBT.  In 

22%  (2/9)  of  the  series,  convergence  required  the  same  number  of  iterations  as 
DRES  PROBT,  but  in  66%  (6/9)  of  the  series,  PROBIT  converged  using  on  the  average 
1/3  of  the  number  of  iterations  of  DRES  PROBT.  In  one  series,  DRES  PROBT  failed 
to  converge  whereas  PROBIT  succeeded. 

4.2  Comparison  of  "PROBIT",  "SAS"  and  "SIPS11 


PROBIT  results  were  tested  against  two  commerical  PROBIT  analysis  programs, 
SAS  and  S103,  using  data  sets  listed  in  Table  5.  The  data  sets  included  several 
series  with  single  and  multiple  null  and  100%  responses.  Such  data,  which  are 
incompatible  with  DRES  PROBT,  provide  a  severe  test  of  program  robustness. 


& 

$ 

8 

w 


Estimates  of  representative  quantities  obtained  for  PROBIT  (Table  6)  were 
very  similar  to  those  obtained  for  S103  and  slightly  different  from  those  of  SAS. 
In  any  case,  the  differences  in  values  form  one  program  to  another  were  not 
statistically  significant.  S103  did  not  have  a  convergence  criterion  but  rather 
ran  for  a  fixed  number  of  iterations,  in  this  case  20,  which  is  far  more  than  is 
normally  required.  This  suggests  that  SAS  may  have  been  one  or  two  iterations 
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away  from  true  convergence.  In  33%  (3/9)  of  the  series,  PROBIT  required  one  less 
iteration  than  SAS  and  in  one  series,  SAS  required  one  less  iteration  than 
PROBIT.  The  biggest  differences  between  all  three  programs  occurred  in  the 
calculation  of  ED99  and  the  fiducial  limits,  being  the  most  pronounced  for  the 
ED99.  This  is  because  fiducial  limits  are  quite  sensitive  to  slight  errors  in 
estimation  of  the  regression  slope,  and  because  any  error  in  calculating  the 
regression  slope  is  compounded  as  one  estimates  quantities  progressively  further 
from  the  mean  of  the  data. 

4.3  "PROBIT"  Transportability 

Finally,  to  test  transportability  of  the  code,  the  PROBIT  source  file  was 
transferred  on  9  track  digital  magnetic  tape  to  the  ODG  Digital  Equipment  VAX 
11/780,  recompiled,  and  relinked.  The  data  sets  of  Table  3  and  5  were  reanalysed 
and  the  results  were  found  to  be  identical  to  those  from  the  Honeywell  computer. 


5.0  DISCUSSION  AND  CONCLUSIONS 


Clearly,  PROBIT  performs  much  better  than  the  old  DRES  PROBT  program  and  is 
comparable  to  the  probit  programs  contained  in  the  two  commerical  packages,  S103 
and  SAS.  It  may,  indeed,  actually  perform  slightly  better  than  SAS. 


Commerical  statistical  package  programs  are  provided  in  an  executable  image 
form,  which  means  that  they  cannot  be  modified.  This  can  be  a  major  drawback, 
since  output  must  then  be  tailored  to  satisfy  a  large  number  of  scientists  of 
differing  disciplines.  The  output  of  the  PROBIT  program  is  currently  configured 
to  the  needs  of  a  much  smaller  number  of  DRES  scientists  with  relatively  similar 
interests.  Further,  it  can  be  modified  as  needs  change. 


Finally,  commerical  statistical  packages  are  exceedingly  expensive. 
Curently,  they  are  rented  for  upwards  of  $1000  per  month  or  are  purchased  for 
between  $10000  and  $25000.  The  development  of  PROBIT  cost  approximately  $8000. 

In  light  of  these  facts,  the  development  of  the  PROBIT  program  has  been  a 
cost  effective  and  worthwhile  endeavour. 
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APPENDIX  1 


Iterative  Estimation  of  Effective  Dose 


A  complete  discussion  of  the  theory  of  probit  analysis  is  beyond  the  scope 
of  this  report  and  may  be  found  in  the  definitive  text  of  Finney  (1971)  and 
Eisenhart  et  al .  (1947).  The  present  technique  is  a  modification  of  those  of  a 
number  of  authors,  primarily  Finney  (1971).  We  present  here  sufficient  theory  to 
make  the  method  clear. 

The  probability  that  an  organism  will  respond,  (e.g.,  die)  if  given  a  dose, 
x‘,  of  a  drug  can  be  represented  by  assuming  that: 


1.  Each  organism  is  characterized  by  a  definite  sensitivity,  s',  (the 
minimum  dose  at  which  it  will  respond), 

2.  sensitivities  are  log-normally  distributed  and 


3.  organism  responses  are  uncorrelated. 

The  justification  and  validity  of  these  assumptions  have  been  discussed 
extensively  (Finney,  1971),  but  for  most  biological  experiments,  the  assumptions 
appear  valid. 


Thus  if  we  define 


s  =  1  n(  s ' ) 


x  =  ln(x')  (x  is  called  a  "dose  metameter") 


then  the  probability  density  function  of  s,  p(s),  is  given  by: 


p(s)  =  o-i  (2*)-i'2  exp  {-(s-p)2/202} 
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The  proportion  of  animals  that  do  not  respond  for  an  Infinite  sample, 
denoted  Q(x),  Is  given  by 


Q(x)  =  o'1  (2n)-i/2  exp  J$=x  exp  {-(s-„) 2/202}  ds  (3) 

The  proportion  that  responds,  called  the  'response',  for  a  dose  metameter  x, 
Is  given  by 


P(x)  =  1  -  Q(x)  (4) 

In  practice,  the  sample  size  is  not  infinite  and  instead  there  is  a  set  of 
trials  (i.e.  treatment  groups)  {i  =  1,  2,  ...n}  in  which  a  given  dose  metameter, 
x. ,  yields  a  number  of  responses,  r. ,  for  n^  organisms.  If  we  can  estimate  p  and 
a  from  the  data,  the  entire  distribution  is  known  and  hence  the  effective  dose 
for  any  percentile  can  be  calculated. 

We  define  =  r^/n^  and  Qi  =  1  -  P^ .  It  can  be  easily  shown  that  P^  and 
are  samples  from  a  binomial  distribution.  Q..  is  an  estimate  of  Q(xi )  where 


Q(xj )  =  a'*  (2*)-i/2  /s=x  exp  {-(s-M)2/202j  ds 


»  (2x)-i/2  /s=Y{x_}  exp  {-s 2/2}  ds 


where  Y(x^)  =  (x^  -  y)/o  (6) 

Y  in  this  form  is  said  to  be  a  'normal  equivalent  deviate',  or  'NED',  A 
probit,  Y',  is  defined  as: 

Y'  =  Y  +  5  (7) 
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can  be  estimated  by  y^  where 


Qi  ■  (2n)-1/2  /  exp  {-s 2/2 }  ds 


(8) 


y^  may  be  obtained  from  by  means  of  tables  or,  more  appropriately  for  computer 
usage,  by  means  of  a  transformation  (y^  =  (Q^ ) )  such  as  is  found  in  Hastings 
(1955).  (It  should  be  noted  at  this  point  that  if  the  natural  response  rate,  PQ, 
defined  as  the  fraction  of  organisms  that  respond  in  the  absence  of  stimulus  for 
an  infinite  population,  is  significant  with  respect  to  the  response  rates  when 
the  stimulus  is  present,  then  it  is  straightforward  to  show  that  the  estimate  of 
the  response  due  to  the  stimulus  for  trial  i,  denoted  P!  is  given  by: 


p;  =  <pi  -  Pc*/*1  -  poJ 
and  Q.  *  1  -  P! 


(9) 

(10) 


Q.j  is  then  used  in  equation  (8)  and  subsequent  equations. 

We  thus  obtain  a  set  of  pairs  of  data  (x^ ,  y^  )  which  should  be  linearly 

related  if  the  assumptions  listed  previously  are  correct,  since: 

Y(x.)  =  Aq  +  AlX.  (11) 

where  Aq  =  -po"1 

and  Aj  =  a"1 

Using  equation  (11)  as  a  model,  the  coefficients  AQ,  At  may  be  estimated  by 
performing  a  weighted  linear  least  squares  regression  of  y^  to  x.. 

Assuming  appropriate  weights,  w^ ,  for  each  trial  i,  are  known  (we  shall 
return  to  this  point),  the  estimate  of  Ap  called  the  'sensitivity',  is  denoted 
a,,  and  it  and  its  estimated  variance,  S3  2,  are  given  by: 

1  a  i 
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a.  =  S  2/S  2 
1  xy  x 


s  2=(xA,2  v-1)/S2 


(12) 

(13) 


The  variance  element  S2  and  covariance  element  S  2  are  given  by: 


xy 


V  *  <1^  -  <1^  niwixi^2  114) 

V  “  (£j  WiV  -  <|"j  "tVt>  1151 


and  xv2  is  an  estimate  of  the  chi-square  parameter  for  the  regression 
(number  of  degrees  of  freedom,  v  =  n-2)  given  by: 


2  = 


V  -  a,S  2 
y  1  xy 


(16) 


where 


y  *  ‘f.!  Vi|(f=1  niv?>  -  wi’2 


(17) 


(xv2v-i  is  called  an  estimate  of  the  "reduced  chi  square" ) 


We  also  define: 


x  =  if  n1wixi )/();"  n.w.) 
i=l  11  i=l  11 

y  *  (f  n.w.y. )/( £n  n.w.) 
i=l  111  i=l  1  1 


(18) 

(19) 


Then  we  can  estimate  Y ( x ^ )  by  y(x^ ) 


yfx^  =  aL  (xi  -  x)  +  y 


(20) 
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The  parent  probability  density  function  from  which  y..  is  a  sample  is  not 
actually  known.  Generally,  it  is  not  a  normal  distribution  and  thus  a  least 
squares  estimate  of  would  not  be  a  maximum  likelihood  estimate.  To  a  first 
approximation,  this  is  equivalent  to  saying  that  the  variance  of  y(x^)  is  a 
function  of  y(x^).  Thus  the  problem  can  be  circumvented  by  weighting  each 
sample  (Bevington,  1969)  but  the  appropriate  weights  to  use,  w. ,  are  not  a  priori 
known.  To  estimate  the  proper  weights,  an  iterative  procedure  (given  for  the  kth 
iteration  where  iteration  number  is  shown  as  a  bracketted  superscript)  is  used): 

1.  calculate  a,'k».  S\™ .  S^».  s/<k>,  :?/<*>,  x  y  <k», 

y(x1-)(k)  using  (12)  -  (20) 

2.  calculate  P.(k)  =  T2(y(x.)(k),  (21) 

i(x.)(k>  ^ 

=  (2k)-W2  /  exp  {-s2/2 }  ds 

S  =  -oo 

3.  calculate  =  1  -  P.(k)  (22) 

4.  calculate  Z.(k)  =  (2n)“i/2  exp  {-y.(k)/2}  (23) 

5.  calculate  Y.(k+1)  =  y(k)(x.)  +  (Q.(k)  -  Q.)/Z.(k)  (24) 

6.  w.(k+1)  =  Z12(k)/{Q.(k)[(l-P0)-l  -  Q.(k)]}  (25) 

7.  go  to  1  for  k+1  iteration. 

To  begin  the  iteration  process,  the  following  are  used  as  initial  values: 

Q.  ^  =  Q. 

yi 

V11*! 
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YiU)  =  Tl(Qi(1)) 

The  iteration  is  repeated  until  convergence  occurs  or  until  a  maximum  number 
of  iterations  have  been  performed.  The  criterion  6  =  -  a^'^l  may  be 

used  to  test  for  convergence. 

Failure  to  converge  is  usually  caused  by  one  of  two  conditions: 

1.  The  data  are  heterogeneous  (see  Appendix  2). 

2.  The  experimental  design  is  not  adequate,  usually  because  the  input  data 
contains  several  trials  with  null  and/or  100%  oberved  responses. 

The  justification  for  the  method  is  not  obvious  and  is  beyond  the  scope  of 
this  report,  but  it  can  be  shown  that  convergence  is  assured  and  that  maximum 
likelihood  estimates  of  the  parameters  are  obtained,  provided  that  the  previously 
discussed  assumptions  are  valid  (Finney,  1971).  Reasons  for  estimating  A0,  Aj 
instead  of  a,  u  are  are  also  discussed. 

When  iteration  is  terminated,  we  can  determine  the  effective  dose  at  the  a 
percentile  level  or  ED ^ ;  that  is  the  dose  at  which  A  percent  of  an  infinite 
population  of  the  organisms  would  respond.  The  ED^,  denoted  symbolically  as  x^, 
is  calculated  as  follows: 


Px  =  A/100 


>  x  =  T.(l-(V 


Xx  *  a;i  (yx  -  y)  +  x 


where  a,,  y,  x,  are  the  values  of  a,  ,  y*k',  x^k*  after  the  last  iteration 


EDX  =  x;  =  exp  (*  ) 
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Error  analysis  must  be  carried  out  after  the  effective  dose  estimation  in 
order  to  determine  the  following: 

1.  the  uncertainties  of  the  estimated  parameters,  and 

2.  the  validity  of  the  assumptions  concerning  the  data  as  discussed  in 
Appendix  1. 

As  in  Appendix  1,  only  an  outline  of  the  error  analysis  will  be  given  and 
the  interested  reader  should  refer  to  any  standard  text  on  probit  analysis  (e.g., 
Finney,  1971)  for  theroretical  justifications  and  a  discussion  of  limitations. 

xa2(v)  is  defined  as  the  value  of  a  variable  from  a  chi -squared  distribution 
such  that 

7  pv  ( x2>  dx2  “  a,  (30) 

x2  =  x2  (  u) 

where  P ^ ( x 2 )  is  the  probability  density  function  of  a  chi-square 
distribution  of  v  degrees  of  freedom.  PJx2)  may  be  calculated  from  an  equation 
or  from  tables  in  any  standard  statistical  text. 

If  x2  >  x2  (  v)  (see  (16)),  it  is  said  that  the  data  is  "heterogenous".  A 
heterogeneity  factor,  h,  is  defined  as 

h  =  l2/v  if  \2  >  x2U>  (3D 

=  1  otherwise. 

The  significance  of  the  heterogeneity  factor  is  frequently  misunderstood. 
Data  will  be  heterogenous  if  the  iterative  regression  analysis  of  Appendix  1  does 
not  correspond  to  a  maximum  likelihood  estimation.  This  may  be  caused  by: 
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1.  correlation  of  response  within  a  trial  due  to  a  lack  of  independence  of 

subjects  within  the  treatment  group  (This  may  lead  to  increased 

dispersion  about  the  regression  line  and  a  large  value  of  x2*  If  this 

is  the  case,  h  can  be  regarded  as  a  factor  by  which  all  weights  have 

been  underestimated  since  for  homogeneous  data  x2  *«•  The  variance  and 

covariances  S2,  S2,  S2  can  then  be  multiplied  by  h.),  or 
x  y  xy 

2.  an  incorrect  mathematical  model.  (One  cannot  compensate  for  such 
deviations  within  the  framework  of  the  analysis.) 

It  is  generally  impossible  to  separate  the  first  type  of  deviation  (random 
deviation)  from  the  second  type  (systematic  deviation). 

The  best  advice  is  to  treat  all  data  which  are  found  to  be  heterogeneous 
with  great  caution,  and  never  to  multiply  variances  by  h  as  a  "cure"  for 
heterogenei ty  unless  one  can  be  sure  that  systematic  deviations  do  not  occur. 
Inspection  of  the  regression  data  and  line  is  essential.  In  the  present  program, 
a  warning  message  is  issued  whenever  heterogeneous  data  occur,  but  the 
multiplication  of  variances  by  h  is  still  carried  out.  Frequently  it  is  found 
that  if  convergence  occurs  the  data  will  be  homogeneous. 

A  final  check  of  goodness  of  fit  can  be  obtained  by  computing  the  expected 
frequence  of  response,  p^  from  the  final  model  after  regression,  for  each  trial. 
If  a  number  of  expected  frequencies  (or  their  complements)  are  very  small,  one 
should  suspect  any  large  value  of  x2* 

Let  t  (v)  be  a  variable  from  a  Student’s  t  distribution  (see  any  standard 
a 

statistics  text)  for  v  degrees  of  freedom  such  that 
t  U) 

ja  P(t)dt  =  1  -  u  (32) 

-t  U) 

a 

where  P(t)  is  the  probability  density  function  of  t.  We  define  a  quantity  g  (the 
g  factor)  as 
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g  =  t2  S  2/a  ?  =  ht  2/  (ajS  2) 
u  a  j  1  a  xy 

(33) 

where  t2  =  t2(v)  h  >  1 

(34) 

a  a 

=  t2(»)  h  =  1. 

a 


Fiducial  limits  for  x  (the  range  of  x  over  which  the  true  x  is  certain  to 

A  A 

lie  at  the  a  confidence  level)  is  given  by: 


f 


A 


g 

X  +  —  (x  -x) 
A  1-g  x 


(l-g)  'Vx)2 


l  Vi  Sx 

1=1  J 


1/2 


a  j  ( l-g) 


(35) 


where  the  plus  sign  corresponds  to  the  upper  limit,  f  ,  and  the  minus  sign 
corresponds  to  the  lower  limit,  fj  . 

The  fiducial  or  Fieller  limits  on  x^,  i.e.,  the  ED^  are  then 


f!  *  exi>  ffxi 


(36) 


It  can  be  readily  shown  that 


g  =  t2  v*1  (r-2  -  1)  (37) 

a 

where  r  is  the  correlation  coefficient  (Bevington,  1969)  for  the  regression. 
Since  -1  <  r  <  1,  then  g  >  0. 


Clearly  if  the  regression  data  are  highly  correlated  or  anti  correlated 
(r  *  ±1),  g  will  be  approximately  zero  and  the  fiducial  limits  are  well  defined. 
If  g  >  1,  ie,  r2  <  (vf2  +  l)-i  then  the  lower  fiducial  limit  is  greater  than  the 
upper,  a  nonsensical  result.  Detailed  analysis  shows  that  when  this  occurs,  the 
range  of  x  over  which  the  true  x  is  certain  to  occur  at  the  a  confidence  level 

A 

is  outside  the  range  given  by  the  two  limits.  Values  of  g  for  "well  behaved" 
data  are  typically  0.1  and  seldom  exceed  0.4  (Finney,  1971).  For  values  of  g 
larger  than  0.4,  the  experimenter  should  view  the  experiment  with  some 
skepticism. 
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Operating  Instructions 


The  probit  analysis  program  PROBIT  is  currently  running  on  the  ORES 
Honeywell  DPS-8  and  the  ODG  Digital  Equipment  VAX  11/780  computers.  Procedures 
for  operating  the  program  on  both  computers  follow. 

Procedure  for  Interactive  Use  on  the  Honeywell  Computer: 

(1)  Log  on  the  terminal  with  the  following  account  number  and  password: 

B 1000200, STATS, DRES.  Press  the  RETURN  key.  Enter  N  to  profile 

change. 

(2)  A  menu  consisting  of  several  statistical  analysis  options  will  appear. 
Find  the  number  corresponding  to  the  probit  analysis  option  (2).  Enter 
the  number  and  press  the  RETURN  key. 

(3)  User  will  be  prompted  with  question: 

Do  you  have  an  input  file  already  created?  (Y/N) 

(4)  If  answer  to  #3  is  N,  then  user  is  prompted  with: 

Do  you  want  to  create  an  input  file  before  running  PROBIT?  (Y/N) 

(5)  If  answer  to  #4  is  N,  then  PROBIT  program  will  prompt  user  to  input 
data  while  the  program  is  running. 

(6)  If  answer  to  #4  is  Y,  then  a  paragraph  of  steps  of  explanation 
detailing  how  to  build  a  file  and  the  format  for  the  data  is  described. 
A  sample  file  is  also  provided  as  an  example. 

(7)  If  answer  to  #3  is  Y,  then  user  is  prompted  with: 

'Enter  file  name' 

After  that  entry,  PROBIT  will  run  successfully. 
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(8)  Perform  Interactive  probit  analysis  for  a  data  set.  The  computer, 

using  simple  instructions,  asks  the  user  for  the  appropriate  input 
data.  Numbers  on  a  line  are  separated  by  a  blank  or  comma.  If 

incorrect  information  has  been  entered  and  the  RETURN  key  has  not  be 
pressed,  the  BACKSPACE  key  may  be  used  to  delete  the  erroneous  data 
which  may  then  be  re-entered.  If  incorrect  information  has  been 

entered  and  the  RETURN  key  has  been  pressed,  the  ESC  Y  key  must  be 
pressed  an  all  previous  data  must  be  re-entered. 

(9)  After  analysis  of  each  data  set,  PROBIT  will  ask  if  you  wish  to  exit 

the  program.  If  'yes'  is  chosen,  the  main  menu  will  be  returned  and 

you  may  carry  out  other  analyses  or  exit  from  the  statistics  package. 

(10)  Afterwards,  the  output  may  be  picked  up  from  the  line  printer. 

Procedure  for  the  VAX  11/780  Computer: 

(1)  Log  on  the  terminal  under  your  account. 

The  terminal  should  previously  be  set  for  132  columns  of  output. 

(2)  If  input  and  output  are  to  be  via  your  terminal  type  the  following: 

@  DISKSMRCGl: [JEM.STATSjPROBITTERM 


(3)  Alternatively,  the  user  can  enter  the  input  from  and  place  the  output 
in  data  files  by  typing: 

@  DISKSMRCGl: [JEM.STATS]PROBITFILE 

It  is  assumed  that  the  data  to  be  entered  have  already  been  placed  in  a 
file  using  the  EDT  editor.  The  input  data  format  is  demonstrated  in 
Appendix  4.  You  will  be  prompted  for  the  file  name  and  for  an  output 
filename.  The  output  file  need  not  initially  exist. 


(4) 


Perform  probit  analysis  for  a  data  set.  When  using  PROBITTERM  the 
interactive  process  is  the  same  as  for  the  Honeywell  computer  except 
that  the  RUBOUT  and  CTRL  Y  keys  are  used  for  data  correction  instead  of 
the  BACKSPACE  and  ESC  Y  keys  respectively. 
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(5)  After  analysis  of  each  data  set,  the  computer  will  ask  the  user  whether 
he/she  wishes  to  terminate  analysis.  If  'yes'  is  chosen,  you  will  be 
returned  to  the  monitor  (DCL)  level. 

(6)  To  logoff,  type  LOG. 


UNCLASSIFIED 


UNCLASSIFIED 


APPENDIX  4 


SAMPLE  INPUT  DATA  FILE  AND 
SAMPLE  OUTPUT  OF  PROGRAM  PROBIT 
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SAMPLE  INPUT  DATA  FILE  : 


SERIES  9  GROUP  2 

4,50 

YES 

10.  585.  10. 

10.  526.5  6. 

10.  473.9  2. 

10.  426.5  3. 

0.95 

0.0  0.001 


SAMPLE  PROBIT  OUTPUT: 
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****  MAINLINE  PROGRAM  **** 

-  A  PROBIT  ANALYSIS  PROGRAM  USED  MAINLY  TO  DETERMINE  THE  - 

PROBABILITY  THAT  AN  ANIMAL  WILL  RESPOND  IF  GIVEN  A  CERTAIN 
EFFECTIVE  DOSAGE  OF  A  DRUG.  ALSO  USED  TO  CALCULATE  THE  REGRESSION 
SLOPE,  PROBITS,  HETEROGENEITY  FACTOR,  AND  FIDUCIAL  LIMITS  FOR  THE 
GIVEN  DATA. 


****  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  *♦** 

-  FUNCTION  XNED(Q)  - 

CALLED  TO  CALCULATE  THE  NED  (NORMAL  EQUIVALENT  DEVIATE)  VALUE 
FOR  A  PROBABILITY  OF  NO  RESPONSE,  Q. 

-  FUNCTION  RESP(Y)  - 

CALLED  TO  CALCULATE  THE  PROBABILITY  OF  RESPONSE,  P,  GIVEN  THE 
NED  VALUE,  Y. 

-  FUNCTION  PCHISQ (CHISQR , NFREE)  - 

CALLED  TO  EVALUATE  THE  PROBABILITY  FOR  EXCEEDING  CHI  SQUARE, 
GIVEN  THE  COMPARISON  VALUE  OF  REDUCED  CHI  SQUARE  (CHISQR)  AND 
THE  NUMBER  OF  DEGREES  OF  FREEDOM  (NFREE) . 

-  FUNCTION  TPROB (FREE , ALPHA)  - 

USED  TO  FIND  THE  VALUE  OF  »T»  FROM  A  T  DISTRIBUTION  TABLE  OF 
'FREE'  DEGREES  OF  FREEDOM,  WHICH  IS  EXCEEDED  WITH  PROBABILITY 
OF  ’ALPHA*  (2  SIDED  TEST). 


****  DESCRIPTION  OF  VARIABLES  **** 

-  INPUT  VARIABLES  - 

ID/  ALPHABETIC;  IDENTIFICATION  TITLE  (PRINTED  ON  EACH  PAGE) 
NT/  INTEGER;  NUMBER  OF  GROUPS 

MAXCNT/  INTEGER;  MAXIMUM  LOOP  COUNT  FOR  ITERATIVE  REGRESSION 
IYNSW/  ALPHABETIC;  SWITCH  FOR  OUTPUT  OF  INTERMEDIATE  RESULTS 
IY  -  YES 
IN  -  NO 

TN/  REAL;  NUMBER  OF  ANIMALS  IN  THE  GROUP 
XIN/  REAL;  DOSE  GIVEN  IN  THE  GROUP 
R/  REAL;  NUMBER  OF  RESPONSES  IN  THE  GROUP 
ALPHA/  REAL;  CONFIDENCE  LEVEL  OF  FIDUCIAL  LIMITS 
DEFAULT  =0.95 

C/  REAL;  FRACTION  OF  NATURAL  MORTALITY 
DEFAULT  =0.0 

DELTA/  REAL;  CONVERGENCE  CRITERION 
DEFAULT  =0.001 

-  CALCULATED  VARIABLES  - 

I, J,K/  INTEGER;  LOOP  COUNTERS 

X/  REAL;  DOSAGES  (NATURAL  LOG  OF  THE  DOSE, XIN) 
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0093  C 

0094  C 

0095  C 

0096  C 

0097  C 

0098  C 

0099  C 

0100  C 

0101  C 

0102  C 

0103  C 

0104  C 

0105  C 

0106  C 

0107 
0108 
0109 
0110 
0111 
0112 
0113 
0114 
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P/  REAL;  PROB.  OF  RESPONSE  AT  A  GIVEN  DOSAGE  LEVEL 
Q/  REAL;  PROB.  OF  NO  RESPONSE  AT  A  GIVEN  DOSAGE  LEVEL 
W/  REAL;  CORRECT  WEIGHTING  COEFF.  CORRESPONDING  TO  Y 
Y/  REAL;  EXPECTED  NEDS  FROM  PROVISIONAL  LINE  (EMPIRICAL) 

PLMBDA/  REAL;  ARRAY  FOR  EFFECTIVE  DOSE  LEVELS 
BOLD/  REAL;  OLD  (PREVIOUS)  CALCULATED  REGRESSION  SLOPE 
SNW , SNWX , SNWY , SNWXX , SNWYY , SNWXY/  REAL;  INTERMEDIATE  SUMS 
XBAR/  REAL;  MEAN  OF  THE  DOSAGES 
YBAR/  REAL;  MEAN  OF  THE  WORKING  NED’S 

SXX/  REAL;  VARIANCE  OF  THE  DOSAGES  ] 

SYY/  REAL;  VARIANCE  OF  THE  WORKING  NED’S 

SXY/  REAL;  COVARIANCE  OF  THE  WORKING  NED’S  AND  DOSAGES 

B/  REAL;  CURRENT  REGRESSION  SLOPE  (REGRESSION  COEFFICIENT) 

GREG/  REAL;  Q  FROM  REGRESSION  —  CORRESPONDING  TO  Y 
2/  REAL;  ORDINATE  CORRESPONDING  TO  Y 
CHISQ/  REAL;  VALUE  OF  CHI  SqUARE 
NFREE/  INTEGER;  NUMBER  OF  DEGREES  OF  FREEDOM 
FREE/  REAL;  NUMBER  OF  DEGREES  OF  FREEDOM 
CHISqR/  REAL;  REDUCED  CHI  SqUARE 
H/  REAL;  HETEROGENEITY  FACTOR 
H  =  1.0;  HOMOGENEOUS 
H  <  1.0  OR  H  >  1.0;  HETEROGENOUS 
T/  REAL;  VALUE  OF  ’T*  DISTRIBUTION 
SICMAB/  REAL;  ERROR  IN  SLOPE 

qLMBDA , YLMBDA , XLMBDA , EDLMB/  REAL;  TOLERANCE  LEVELS  OF  VALUES 
C,FACT,FACT1,FACT2,DIV,EX/  REAL;  INTERMEDIATE  VALUES 
FUX#FLX/  REAL;  UPPER  AND  LOWER  FIDUCIAL  LIMITS  FOR  X  (ERRORS) 

FUED.FLED/  REAL;  UPPER  AND  LOWER  FIDUCIAL  LIMITS  FOR  ED 
YR/  REAL;  MAXIMUM  LIKELIHOOD  EXPECTED  PROBITS  (NEDS  -  5) 

PLB/  REAL;  ED  LEVEL  (DOSE  LEVEL  *  100) 

APA/  REAL;  PERCENT  CONFIDENCE  LEVEL 

RP/  REAL;  EMPIRICAL  PERCENTAGE  RESPONSE 

FST/  REAL;  FIRST  COEFFICIENT  OF  REGRESSION  LINE  EqUATION 

YNEW/  REAL;  WORKING  NEDS  (NORMAL  EqUIVALENT  DEVIATES) 

YK/  REAL;  MAXIMUM  LIKELIHOOD  EXPECTED  NEDS 
SIGY/  REAL;  STANDARD  ERROR 

FLY/  REAL;  LOWER  FIDUCIAL  LIMITS  OF  EXPECTED  NUMBER  REACTING 
FJL/  REAL;  UPPER  FIDUCIAL  LIMITS  OF  EXPECTED  NUMBER  REACTING 
RS/  REAL;  MEAN  OF  THE  EXPECTED  NUMBER  REACTING 
DUMMY/  REAL;  ARRAY  FOR  DUMMY  VARIABLE 
SWXY,SWX,SWY/  REAL;  INTERMEDIATE  VALUES 


****  PROGRAMMER  /  DAWN  RIGBY  *♦** 
DATE  /  22/06/82  **** 

****  Updated  /  John  McFee  **** 
****  Date  /  26/02/86  ***+ 


IMPLICIT  REAL*8  (A-H,  0-Z) 

DIMENSION  TN(IOO) ,XIN(100) ,R(100) ,ID(50) ,P(100) ,q(100) ,W(100) 
DIMENSION  X (100) , Y (100) , qREG (100) , Z (100) , YR (100) , YNEW(IOO) 
DIMENSION  qLMBDA (10) ,YLMBDA(10) ,XLMBDA(10) ,PLMBDA(10) ,EDLMB(10) 
DIMENSION  FACT (10) , FACT1 (10) ,FACT2(10) ,FLX(10) ,FUX(10) , FLED (10) 
DIMENSION  PLB (10) ,RP(100) ,YK(100) , FLY (100) ,FUL(100) , SIGY (100) 
DIMENSION  DUMMY (100) ,RS(100) ,FUED(10) 

DATA  IY/’YES  ’/,IN/’NO  ’/ 

UNCLASSIFIED 
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data  i ndev/2/, ioutdev/3/, infodev/6/ 


—  INPUT  DATA  - 


wr i te ( i nf odev , 404) 

FORMAT (IX, ’IF  AN  INCORRECT  VALUE  IS  "ENTERED*  AS  INPUT,  INTERRUPT 
1THE  SYSTEM’, /X, 'IMMEDIATELY  BY  PRESSING  "ESC  Y"  TWICE.’, /X, 

2’BEGIN  THE  INPUT  PROCEDURE  FROM  THE  BEGINNING.’) 
wr i te( i nfodev , 1) 

FORMAT  (IX, ’TYPE  A  SUITABLE  EXPERIMENT  TITLE ’/IX, ’EXCLUDING  DIRECT 

1  MENTION  OF  PROBIT  ANALYSIS.’) 
read(indev,2)ID 

FORMAT (50A1) 
write (inf odev, 3) 

FORMAT  (IX,’ INPUT  THE  NUMBER  OF  GROUPS  AND  THE  MAXIMUM  LOOP  COUNT’ 
1/1X, ’ (IF  UNKNOWN,  TYPE  50  (DEFAULT  =  50)) ; ’/IX, ’SEPARATE  THESE  TWO 

2  INTEGERS  BY  A  BLANK  OR  A  COMMA.’) 
read (i ndev , *)  NT.MAXCNT 

write (inf odev, 4) 

F0RMAT(1X, ’TYPE  EITHER  "YES"  OR  "NO",  DEPENDING  UPON  WHETHER  YOU 
1WANT  THE ’/IX, 'INTERMEDIATE  RESULTS  PRINTED  AS  OUTPUT  OR  NOT.’) 
read ( i ndev , 5) IYNSW 
FORMAT (A4) 
write (inf odev, 6) 

FORMAT (IX, ’INPUT  THE  NUMBER  OF  ANIMALS,  THE  DOSE,  AND  THE  NUMBER 
10F ’/IX, ’RESPONSES  FOR  EACH  GROUP,  ONE  TRIAL  PER  LINE. ’/IX, ’THESE 
2MUST  BE  "REAL"  NUMBERS,  SEPARATED  BY  A  BLANK  OR  A  COMMA.’) 
read ( i nde v , *) (TN (I) , XIN (I) , R (I) ,  1=1 , NT) 
write (inf odev, 7) 

FORMAT (IX, ’INPUT  THE  CONFIDENCE  LEVEL  OF  THE  FIELLER  LIMITS’/1X,  ’ 
1 (DEFAULT=. 95;  REFER  TO  THE  ABOVE  "MAX.  LOOP  COUNT"  EXPLANATION’ ,/X 
2, ’DEFAULT) .  THIS  MUST  BE  A  "REAL"  NUMBER. ’) 
read ( i ndev , *)  ALPHA 
wr i te(i nfodev ,8) 

FORMAT (IX, ’INPUT  THE  FRACTION  OF  NATURAL  MORTALITY  (DEFAULT=0.) ’/ 
11X, ’AND  THE  CONVERGENCE  CRITERION  (DEFAULT=. 001) . ’/IX, ’THESE  MUST 
2B0TH  BE  "REAL"  NUMBERS,  SEPARATED  BY  A  BLANK  OR  A  COMMA.’) 
read (i ndev,*)  C, DELTA 


-  OUTPUT  OF  THE  INITIAL  DATA  - 


write(ioutdev,9  )ID 

FORMAT(’l ’ ,60X, ’PROBIT  ANALYSIS ’/61X,50A1) 
write(i outdev, 10) 

FORMAT (IX ,//,21X, ’INITIAL  INPUT  DATA:’) 
wr i te(i outdev, 11) 

FORMAT (IX,  //,15X, ’GROUP  NUMBER ’, 9X ,’ NUMBER  OF  ANIMALS’, 9X, 
1 ’ DOSE \9X, ’NUMBER  OF  RESPONSES’,/) 

IF (NT. LE. 50) GO  TO  12 
DO  13  1=1,50 

write(ioutdev,14)I,TN(I),XIN(I),R(I) 

CONTINUE 

wr i te( i outdev ,9  )ID 
write(i outdev, 11) 

DO  15  1=51, NT 

wr i te ( i outdev , 14) I , TN (I) ,XIN(I),R(I) 

CONTINUE 

UNCLASSIFIED 
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0172 

GO  TO  16 

0173 

12 

DO  17  1=1, NT 

0174 

write(ioutdev,14)I,TN(I),XIN(I),R(I) 

0175 

14 

FORMAT (18X , 14 , 19X , F7 . 2 , 12X , G9 . 4E2 , 15X , F7 . 2) 

0176 

17 

CONTINUE 

0177 

16 

wr i te ( i outdev , 9) ID 

0178 

wr i te ( i outdev , 18)  NT , MAXCNT , IYNSW 

0179 

18 

F0RMAT(1X,//,12X, ’NUMBER  OF  GROUPS: ’ ,14,//, 12X, ’MAXIMUM  LOOP  COUNT 

0180 

1: ’,14,//, 12X, ’OUTPUT  OF  INTERMEDIATE  RESULTS:  ’,A4) 

0181 

write(i outdev, 19)  ALPHA, C, DELTA 

0182 

19 

F0RMAT(1X,/,12X, ’CONFIDENCE  LEVEL  OF  FIELLER  LIMITS: * ,F8. 4,//, 

0183 

112X, ’FRACTION  OF  NATURAL  MORTALITY: ’,F8.5//12X, ’CONVERGENCE  CRITER 

0184 

2I0N: ’ ,F8 .5) 

0185 

C 

0186 

C 

-  INITIAL  DATA  VALUES  - 

0187 

C 

ASSUME  SENSITIVITIES  ARE  LOG  NORMALLY  DISTRIBUTED 

0188 

C 

0189 

write(i outdev, 20  ) 

0190 

20 

F0RMAT(1X//////21X, ’INITIAL  CALCULATED  VALUES:’) 

0191 

write(i outdev, 21  ) 

0192 

21 

FORMAT (1X//10X, ’GROUP  NUMBER’, 9X, ’DOSAGES’ ,9X, ’EMPIRICAL  RESPONSES 
1 ’ ,9X, ’EMPIRICAL  NEDS’, 9X, ’WEIGHTS’,/) 

0193 

0194 

DO  22  1=1, NT 

0195 

X  (I) =DLOG (XIN (I) ) 

0196 

P(I)=R(I)/TN(I) 

0197 

P(I)  =  (P(I)-C)/(1.0-C) 

0198 

Q  (I)  =i  -  o— p  (i) 

0199 

W(I)=1.0 

0200 

IF ( (P (I) .EQ.O.O) .or . (q(I) .eq.O.O))  W(I)=0.0 

0201 

Y  (I)  =XNED  (Q.(I) ) 

0202 

write(i outdev, 23  )I,X(I)  ,P(I)  ,Y(I)  ,W(I) 

0203 

23 

FORMAT (14X,I4,12X,E12.5,10X,G9.4E2,16X,E12.5,11X,F7.3) 

0204 

22 

CONTINUE 

0205 

BOLD=1.D*30 

0206 

KLMBDA  =  5 

0207 

PLMBDA(1)=0 .01 

0208 

PLMBDA(2)=0.50 

0209 

PLMBDA (3) =0 . 90 

0210 

PLMBDA(4)=0 .95 

0211 

PLMBDA (5) =0.99 

0212 

C 

0213 

C 

-  REGRESSION  LOOP  - 

0214 

C 

0215 

C 

—  FORM  SUMS  — 

0216 

C 

0217 

IF (IYNSW. EQ. IN)  GO  TO  24 

0218 

wr i te ( i outdev , 25) ID 

FORMAT (’l’,60X, ’PROBIT  ANALYSIS ’/61X,50A1) 

0219 

25 

0220 

wr i te ( i outdev , 26) 

0221 

26 

FORMAT (1X,//,21X, 'INTERMEDIATE  RESULTS: ’) 

0222 

24 

DO  27  J=l, MAXCNT 

0223 

SNW=0 . 

0224 

SNWX=0. 

0225 

SNWY=0. 

0226 

SNWXX=0 . 

0227 

SNWYY=0. 

0228 

SNWXY=0 . 

UNCLASSIFIED 
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m 
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0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

28 

0237 

C 

0238 

C 

0239 

C 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

30 

0247 

31 

0248 

0249 

29 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 

0261 

0262 

0263 

34 

0264 

33 

0265 

0266 

0267 

35 

0268 

0269 

0270 

37 

0271 

0272 

39 

0273 

40 

0274 

0275 

38 

0276 

42 

0277 

27 

0278 

43 

0279 

C 

0280 

C 

0281 

c 

0282 

36 

0283 

44 

0284 

0285 

45 

UNCLASSIFIED 
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DO  28  1=1, NT 
SNW=SNW  «TN(I)  *  W(I) 

SNWX=SNWX  ♦TN(I)  *  W(I)  *  X(I) 

SNWY=SNWY  ♦TN(I)  *  W(I)  *  Y(I) 
SNWXX=SNWXX  ♦TN(I)  *  W(I)  *  X(I)**2 
SNWYY=SNWYY  ♦TN(I)  *  W(I)  *  Y(I)**2 
SNWXY=SNWXY  *TN(I)  *  W(I)  *  X(I)  *  Y(I) 
CONTINUE 


—  CALCULATE  STATISTICS  OF  DATA  — 


XBAR=SNWX/SNW 
YBAR=SNWY/SNW 
SXX=SNWXX/SNW  -  XBAR**2 
SYY=SNWYY/SNW  -  YBAR**2 
SXY=SNWXY/SNW  -  XBAR*YBAR 
IF (SXX-0 . 0) 29 , 30 , 29 
write(ioutdev,31) 

FORMAT (’1’, IX, ’DOSES  HAVE  ZERO  VARIANCE— REGRESSION  IMPOSSIBLE!’) 

GO  TO  32 

B=SXY/SXX 

CHISQ=  (SYY  -  B*SXY)  ♦  SNW 
DO  33  1=1, NT 

Y (I)=B  *  (X(I)  -  XBAR)  ♦  YBAR 
IF (Y (I)  .GT.  5.0)  Y (I)  =  5.0 
IF(Y(I).LT.  -5.0)  Y (I)  =  -5.0 
QREG(I)  =  1.0  -  RESP (Y (I)) 

EX=(-(Y(I)**2)/2.) 

Z (I) = (DEXP (EX) ) /DSQRT (2 . *3 . 1415926535D0) 

YNEW (I) =Y (I) ♦ (QREG (I) -  Q(I))/Z(I) 

DIV  =  (QREG (I)  *  ((1-/(1-  O)  -  QREG (I))) 

IF (DABS  (DIV) . LT . 1 . 00D-10) GO  TO  34 
W(I)  =  (Z(I)**2)/(DIV) 

GO  TO  33 
W(I)=0.0 
CONTINUE 

IF(IYNSW.EQ.IN)  GO  TO  35 
GO  TO  36 
DO  37  1=1, NT 

YK(I)  =  Y(I) 

Y  (I)  =  YNEW(I) 

CONTINUE 

IF  (  J-MAXCNT) 38,38, 39 
wr i te( ioutdev ,40) 

FORMAT (’1’, IX, ’WARNING  —  REGRESSION  MAY  NOT  BE  CONVERGENT’) 

GO  TO  41 

IF (DELTA  -  DABS (B-BOLD) ) 42 , 43, 43 
BOLD  =  B 
CONTINUE 
GO  TO  41 


-  OPTIONAL  OUTPUT  OF  INTERMEDIATE  RESULTS 


wr i te ( i outdev , 44) J 

FORMAT (IX, ///,2X, ’LOOP  NUMBER: ’,14) 

write(ioutdev,45) 

FORMAT (1X//4X, ’REGRESSION’, 6X, ’VARIANCE  0F’,6X,’VAR.  OF  PREVIOUS’, 
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0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 

0307 

0308 

0309 

0310 

0311 

0312 

0313 

0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 

0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 
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16X/C0VAR.  OF  PREVIOUS ’,6X, ’MEAN  OF’, 6X, ’MEAN  OF  PREVIOUS’ ,6X, 
2*CHISQ’/4X, ’COEFF. (SLOPE) ',3X, ’DOSAGES ’, 10X, ’LOOP  WORKING  NEDS’, 
35X, 'LOOP  WORKING  NEDS’ ,7X, ’DOSAGES' ,6X, ’LOOP  WORKING  NEDS’/4X, 

4 ’  — B’ ,  13X,  ’— SXX ’ ,  12X,  ’ — SYY ’ ,  17X,  ’AND  DOSAGES— SXY ’ , 8X ,  ’ - -XBAR ’ , 
57X, ’--YBAR’ ,/) 

wr i te ( i outde v , 46) B , SXX , SYY , SXY, XBAR , YBAR, CHISQ 
FORMAT (4X , E12 . 5,4X, E12 . 5,5X,E12 . 5, 12X,E12 .5, 7X, E12 .5,5X,E12.5,6X, 
1E12.5) 

wr ite(ioutdev,47  ) 

FORMAT (IX/ //4X , ’ GROUP ’ , 1 IX , ’ DOSAGES ’ , 10X , ’ EXPECTED  NEDS ’ , 9X , ’ WORKI 
1NG’,17X, ’Q  FROM’ ,7X, ’ORDINATE’ ,9X, ’WEIGHT’/4X, ’NUMBER ’, 10X, ’— X (I) 
2 ’ , 1 IX , ’ FROM  PROVISIONAL ’,6X, ’NEDS’, 18X, ’ REGRESSION ’ , 5X , ’ CORRESPOND 
3ING ’ , 4X , ’ CORRESPONDING ’/4X , ’ — I ’ , 30X , ’LINE— Y (I) ’ , 12X, ’ — YNEW(I) ’ , 
413X,  ’--QREG(I)  ’  ,6X,  ’TO  Y(I) —  Z(I)  ’  ,4X,  ’TO  Y(I) — W(I)’/) 

DO  48  1=1, NT 

write(ioutdev,49  ) I , X (I) , Y (I) ,YNEW(I) , QREG (I) , 2 (I) ,W(I) 

FORMAT (4X , 14 , 12X , El 2 . 5 , 5X , E12 . 5 , 12X , E12 . 5 , 7X , E12 . 5 , 5X, E12 . 5 , 6X , 
1E12.5) 

CONTINUE 
GO  TO  35 

-  FURTHER  CALCULATIONS  - 

NFREE  =  NT  -  2 
FREE  =  NFREE 
CHISQR  =  CHISq/FREE 

IF (PCHISq(CHISqR, NFREE). GE.(1. -ALPHA)) GO  TO  51 
H=CHISQR 

T  =  TPROB (FREE, ALPHA) 

GO  TO  52 
H  =  1.0 

T  =  TPROB (900., ALPHA) 

SIGMAB  =DSQRT (H/ (SXX*SNW) ) 

G  =(H  *  T**2)/(B**2  *  SXX  *  SNW) 

DO  105  1=1, NT 

SIGY(I)  =  DSQRT ( (H/SNW)  ♦  (H*(X(I)-XBAR)**2/(SXX  *  SNW))) 

FUL(I)  =  YK(I)  ♦  T  *  SIGY(I) 

FLY(I)  =  YK (I)  -  T  *  SIGY(I) 

DUMMY (I)  =0.0 
CONTINUE 

DO  53  K=1,KLMBDA 
QLMBDA (K) =1 . -  PLMBDA(K) 

YLMBDA (K) =XNED (QLMBDA (K) ) 

XLMBDA (K) =  ( (YLMBDA (K) -YBAR) /B) ♦XBAR 
EDLMB (K) =DEXP (XLMBDA (K) ) 

IF(G.GT.l.O)  GO  TO  53 

-  FIELLER  LIMITS  (ERRORS  IN  ED)  - 

FACT (K) = ( (H* ( 1 . -C) ) /SNW) ♦ (H* (XLMBDA (K) -XBAR) **2) / (SXX*SNW) 

FACTl (K) =DSQRT (FACT (K) ) *  (T/ (B* (1 . -G) ) ) 

FACT2  (K)  =XLMBD A  (K)  ♦  (G/  (1 .  -G) )  *  (XLMBDA (K)  -XBAR) 

FLX(K)=FACT2(K) -FACTl (K) 

FUX (K) =FACT2 (K) ♦FACTl (K) 

FLED (K) =DEXP (FLX (K) ) 

FUED (K) =DEXP (FUX (K) ) 

CONTINUE  UNCLASSIFIED 
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0343  APA  =  ALPHA  *  100. 

0344  C 

0345  C  -  OUTPUT  OF  FINAL  RESULTS  - 

0346  C 

0347  write (ioutdev, 55) ID 

0348  55  FORMAT (’l',60X, ’PROBIT  ANALYSIS ’/61X,50Al) 

0349  write (ioutdev, 56) 

0350  56  F0RMAT(1X,//  ,21X, ’FINAL  OUTPUT  DATA: ’) 

0351  IF (MAXCNT - J) 57 , 58 , 58 

0352  57  wr i te( ioutdev ,59) MAXCNT 

0353  59  FORMAT (IX,  //  ,20X,  ’  !!  WARNING  —  ITERATIVE  REGRESSION  FAILED  TO 

0354  1  CONVERGE  AFTER  ’,13,’  ITERATIONS ’/22X, RESULTS  FOLLOWING  MAY  BE 

0355  2  SUSPECT  ! !  ’) 

0356  GO  TO  60 

0357  58  wr ite(ioutdev,61) J-l 

0358  61  FORMAT (IX,  //  ,20X,  ’  **  ITERATIVE  REGRESSION  CONVERGED  AFTER  ’,13 

0359  1,’  ITERATIONS  **  ’) 

0360  60  write (ioutdev, 62) 

0361  62  FORMAT (IX , ////, 6X ,’ GROUP ’, 6X ,’ DOSE ’, 6X , ’ DOSAGES ’ , 6X , ’ EMPIRICAL ’ , 6X 

0362  1, 'MAXIMUM  LIKELIHOOD ’ ,6X, ’MAXIMUM  LIKELIHOOD’ ,6X, ’STANDARD’ ,6X, 

0363  2 ’WEIGHT’) 

0364  write (ioutdev, 63) 

0365  63  FORMAT (6X, ’NUMBER \15X, ’ (LN  0F*,7X,’%  RESPONSE’ ,5X, ’EXPECTED  PROBI 

0366  1TS’,8X, ’EXPECTED  NEDS’ , 11X, ’ERR0R’/27X, ’DOSE) ’//) 

0367  IF (NT. LE. 50) GO  TO  64 

0368  DO  65  1=1,50 

0369  YR(I)=YK(I)  ♦  5. 

0370  RP(I)=P(I) *100.0 

0371  wr i te ( i outdev ,67)1, XIN (I) , X (I) , RP (I) , YR(I) , YK (I) , SIGY (I) , W(I) 

0372  65  CONTINUE 

0373  write (ioutdev, 99) APA 

0374  99  F0RMAT(1X/////,6X, ’GROUP’, 6X, ’DOSE’, 6X, ’DOSAGES’, 6X, ’NUMBER  TESTED 

0375  1’  ,6X,  ’ . - . - . NUMBER  REACTING - 

0376  2 - ’/6X,  ’ NUMBER ’,52X,  ’OBSERVED’ ,6X,  ’ . EXPECTED— 

0377  3- . - . -’/80X,  ’MEAN’  ,6X,  ’  —  ’  ,F6.2,  ’%  FIDUCIAL  LIMITS - 

0378  4—  ’/96X,  ’LOWER ’ ,  11X,  ’UPPER’//) 

0379  DO  101  1=1,50 

0380  RS (I) =RESP (YK (I) )  *  TN(I) 

0381  FLY (I) =RESP (FLY (I) )  *  TN(I) 

0382  FUL (I) =RESP (FUL (I) )  *  TN(I) 

0383  write (ioutdev, 100) I, XIN (I) , X (I) ,TN(I) ,R(I) ,RS(I) ,FLY(I) ,FUL(I) 

0384  100  FORMAT (7X, I3,5X,G9 . 4E2, 2X, E12 .5,6X, F7 .2, 14X,F7 .2,8X,F7 .2,6X,E12.5, 

0385  15X.E12.5) 

0386  101  CONTINUE 

0387  wr i te( i outdev , 55) ID 

0388  write (ioutdev, 62) 

0389  write (ioutdev, 63) 

0390  DO  66  1=51, NT 

0391  YR(I)=YK(I)  ♦  5 . 

0392  RP(I)=P(I) *100.0 

0393  write (ioutdev, 67) I, XIN (I) , X (I) , RP (I) , YR (I) ,YK(I) ,SIGY(I) ,W(I) 

0394  66  CONTINUE 

0395  write (ioutdev, 99) APA 

0396  DO  102  1=51, NT 

0397  RS  (I)  =  RESP(YK(I))  *  TN(I) 

0398  FLY (I) =RESP (FLY (I) )  *  TN(I) 

0399  FUL(I)=RESP(FUL(I) )  *  TN(I) 
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wr i te( i outdev , 100) I,XIN(I) ,X(I) ,TN(I) ,  R(I) ,RS(I) ,FLY (I) ,  FUL(I) 

CONTINUE 

GO  TO  68 

DO  69  1=1, NT 

YR  (I)  =YK  (I)  +5 . 

RP(I)=P(I)  *100 .0 

write(i outdev, 67)1, XIN(I)  ,X(I)  ,RP(I)  ,YR(I)  , YK(I)  ,SIGY(I)  , W(I) 
FORMAT (7X , 13 , 5X , G9 . 4E2 , 2X , E12 . 5 , 3X , F7 . 2 , 10X , E12 . 5 , 12X , E12 . 5 , 7X , 
1E12.5,3X,E12.5) 

CONTINUE 

wr i te ( i outdev , 99) APA 

DO  103  1=1, NT 

RS (I) =RESP (YK (I) )  *  TN(I) 

FLY (I) =RESP (FLY (I) )  *  TN(I) 

FUL (I) =RESP (FUL (I) )  *  TN(I) 

write(i outdev, 100)1, XIN(I) ,X(I) ,TN(I) ,R(I) ,RS(I) ,FLY(I) ,FUL(I) 
CONTINUE 

wr i te ( i outdev , 70) ID 

F0RMAT(’l’,6OX,  ’PROBIT  ANALYSIS 761X,50A1) 

FST=YBAR  -  8  *  XBAR 
write(i outdev, 71  )FST,B 

FORMAT (1X,///,22X, ’THE  EQUATION  FOR  THE  REGRESSION  LINE  IN  NEDS  IS 
1:’/31X,’Y  =  ’ ,G14 .7E2, ’  +’,G12.5E2,’  X’) 
wr i te ( i outdev , 72) B , SIGMAB 

FORMAT (IX,  /  ,22X, ’REGRESSION  SLOPE  (SENSITIVITY)  =  \G12.5E2, 

1’  +0R-  ’,  G12.5E2, ’  (STANDARD  ERROR)*) 
wr i te ( i outdev , 73) CHISQ, NFREE 

FORMAT (1X,/,22X, ’CHI  SQUARE  =’,G12.5E2,’  FOR*, 14,’  DEGREES  OF  FR 
1EED0M’) 

IF(PCHISQ(CHISQR, NFREE). GE. (1 . -ALPHA)) GO  TO  74 
H=CHISQ/FREE 
GO  TO  75 
H=1 .0 

IF(H.EQ.1.0)GO  TO  76 
wr ite(i outdev, 77) H,G,T, NFREE 

FORMAT (IX, /,22X ,’!!!!!  WARNING  —  THE  DATA  IS  HETEROGENOUS:  THE 
1HETER0GENEI I Y  FACTOR  IS  H  =  *,E12.5,’  !!!!!* /22X ,*!!!!!  G  =’,E12.5 
2,’;  T  =’ ,E12.5, *  WITH’, 14,’  DEGREES  OF  FREEDOM.  ! M ! ! ’/22X, ’!!!!! 
3  "TREAT  PROBIT  RESULTS  WITH  EXTREME  CAUTION"  !!!!!’) 

GO  TO  78 

wr i te( i outdev , 79)G,T 

FORMAT (IX, /,22X, ’THE  DATA  IS  HOMOGENEOUS;  THE  HETEROGENEITY  FACTOR 

1  IS  H  =  1.0;  G  =’ ,E12.5,/31X, ’T  =*,E12.5,’  WITH  INFINITE  DEGREES 
20F  FREEDOM’) 

IF(G.LE.l.O)  GO  TO  300 
wr i te ( i outdev , 301) 

FORMAT (//22X, ’!!! i !  WARNING  —  G  >  1.0;  THE  FIELLER  LIMITS  ARE  OF 
1N0  USE  AND  HENCE  NOT  PRINTED  HIM  ’/22X,  ’!!!!!  SEE  FINNEY,  PAGE  79 

2  !!!!!’) 

wr i te ( i outdev , 80) APA 

FORMAT(/////UX, ’EXACT  FIELLER  LIMITS  BELOW  GIVE  THE  \F6.2,’»  FID 
1UCIAL  LIMITS’) 
wr  i te( i outdev , 81) 

FORMAT (//11X,  ’ . . 

1 . ’) 

write(ioutdev,82) 

FORMAT  (/ / 1 3X ,’ ED  ’,  10X ,’  EFFECTIVE  ’,  8X ,  ’ . FIELLER  LIMITS . ’  / 
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0492 

400 
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111X, ’LEVEL  K’,7X, ’DOSE’ , 16X, ’LOWER’, 10X, ’UPPER*) 
wr i te ( i outdev , 83) 

FORMAT (//11X, ’RESPONSE', 6X,’ - IN  NATURAL  LOGARITHMIC  SCALE - 

1-’) 

DO  84  K  =1,KLMBDA 
PLB(K)  =  PLMBDA(K)  *  100. 

IF(C.LE.l.O)  GO  TO  302 

wr i te ( i outdev , 303)  PLB (K) , XLMBDA (K) 

FORMAT (12X,  F7 . 2, 8X ,  F8 . 5) 

GO  TO  84 

wr i te ( i outdev , 85)  PLB (K) , XLMBDA (K) , FLX (K) , FUX (K) 

FORMAT (12X , F7 . 2 , 8X , F8 . 5 , 9X , E12 . 5 , 5X , E12 . 5) 

CONTINUE 

write(ioutdev,86) 

FORMAT (/11X, ’RESPONSE ’,6X,’~ . IN  ORIGINAL  SCALE . . 

1-’) 

DO  87  K=1 .KLMBDA 
PLB(K)  =  PLMBDA(K)  ♦  100. 

IF(G.LE.l.O)  GO  TO  304 

wr i te ( i outdev , 305)  PLB (K) , EDLMB (K) 

FORMAT (12X,F7 .2,8X, F8 .3) 

GO  TO  87 

wr i te ( i outdev , 88)  PLB (K) , EDLMB (K) , FLED (K) , FUED (K) 

FORMAT (12X,F7.2,8X,E12.5,9X,E12.5,5X,E12.5) 

CONTINUE 

write(i outdev ,81) 

SWXY  =  SNW  *  SXY 
SWX  =  SNW  *  SXX 
SWY  =  SNW  *  SYY 
write(i outdev ,389) 

F0RMAT(/////11X, ’*♦***  PARAMETERS  NECESSARY  FOR  RELATIVE  POTENCY  T 
1EST  PROGRAM  ***** ’///13X, ’CHISQ’ ,9X, ’SNW’,10X, ’XBAR’,7X, ’SXX*SNW’, 
26X, ’SXY*SNW’,8X, ’DFREE’ ,8X, ’YBAR’,7X, ’SYY*SNW’) 
wr i te ( i outdev , 400)  CHISQ , SNW , XBAR , SWX , SWXY , NFREE , YBAR , SWY 
FORMAT (/10X.E11 .5,3X,F10.5,3X,F10.5,3X,F10.5,3X,F10.5,7X,I4,5X, 
1E11 ,5,3X,Ell .5) 

STOP 

END 
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****  FUNCTION  XNED  **** 

-  THIS  FUNCTION  CALCULATES  THE  NORMAL  EQUIVALENT  DEVIATE  - 

(NED)  VALUE  FOR  A  PROBABILITY  OF  NO  RESPONSE,  Q.  IT  ASSUMES 
THAT  NED(Q)  HAS  A  MAXIMUM  ERROR  OF  ♦  OR  -.0004  OVER  THE  GIVEN 
RANGE  (  MOST  LIKELY,  0  <  Q  <  OR  =  .5  ).  THE  TRANSFORM  USED 
IS  FROM  "APPROXIMATIONS  FOR  DIGITAL  COMPUTERS"  BY  HASTINGS,  P 
RINCETON  UNIVERSITY  PRESS,  PRINCETON,  NJ.,  1955. 


♦***  USAGE  **** 

-  RESULT  =  XNED (Q)  - 


****  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  **♦* 
-  NONE  - 


****  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES  **** 

-  PARAMETERS  PASSED  FROM  THE  MAINLINE  - 

q/  REAL;  PROBABILITY  OF  NO  RESPONSE  (  OF  SURVIVAL) 

-  VARIABLES  - 

ETA/  REAL;  INTERMEDIATE  VALUE 
A0,A1,A2,B1,B2,B3/  REAL;  CONSTANTS 


****  PROGRAMMER  /  DAWN  RIGBY  ***♦ 
****  DATE  /  22/06/82  **** 


FUNCTION  XNED(Q) 

IMPLICIT  REAL*8  (A-H,  0-Z) 
IF(Q.NE.1.0)GO  TO  18 


-  IF  Q  IS  EQUAL  TO  1.0  - 


XNED  =  -5.0 
GO  TO  14 
IF (Q- .5) 2, 2, 4 
IF(Q-0.0)6,7,8 

-  IF  Q  IS  LESS  THAN  OR  EQUAL  TO  0.0  - 


write(ioutdev,10) 

FORMAT (’ 1 ' , IX,  ’SUBROUTINE  XNED— ’/IX, ’NED  VALUE  IS  NOT  DEFINED; 
1Q<0.0’) 

GO  TO  13 
XNED  =5.0 
GO  TO  14 
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0058 

C 

0059 

C 

-  IF  Q  IS  GREATER  THAN  0.5  - 

0060 

C 

0061 

4 

ETA  =DSQRT (-2. *  DLOG(1.0  -  Q)) 

0062 

GO  TO  12 

0063 

C 

0064 

C 

-  IF  q  IS  BETWEEN  0.0  AND  0.5  — 

— 

0065 

C 

0066 

8 

ETA  =DSQRT(-2.*  DLQG(Q)) 

0067 

12 

AO  =  2.515517 

0068 

Al  =  0.802853 

0069 

A2  =  0.010328 

0070 

B1  =  1.432788 

0071 

B2  =  0.189269 

0072 

B3  =  0.001308 

0073 

XNED=  ETA- ( (A0+Al*ETA*A2*ETA**2) / (1 . 

♦Bl*ETA-*-B2*ETA*+2+B3*ETA**3)) 

0074 

IF(q  -  .5)14,14,16 

0075 

C 

0076 

C 

-  IF  Q  IS  GREATER  THAN  0.5  - 

0077 

C 

0078 

16 

XNED  =  -XNED 

0079 

GO  TO  14 

0080 

13 

STOP 

0081 

14 

RETURN 

0082 

END 

UNCLASSIFIED 
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****  FUNCTION  RESP  **** 

-  THIS  FUNCTION  IS  USED  TO  CALCULATE  THE  PROBABILITY  OF  - 

RESPONSE,  P,  GIVEN  THE  NED  VALUE,  Y.  RESP(Y)  HAS  A  MAXIMUM 
ERROR  OF  +  OR  -  3.D-7  AND  USES  THE  TRANSFORMATION  FROM 
"APPROXIMATIONS  FOR  DIGITAL  COMPUTERS"  BY  C.  HASTINGS, 
PRINCETON  UNIVERSITY  PRESS,  PRINCETON,  NJ.,  1955. 


****  USAGE  **** 

-  RESULT  =  RESP(Y)  - 


****  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  ♦*** 
-  NONE  - 


****  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES  **** 

-  PARAMETERS  PASSED  FROM  THE  MAINLINE  - 

Y/  REAL;  NORMAL  EQUIVALENT  DEVIATE  (NED)  VALUE 

-  VARIABLES  - 

X,RES/  REAL;  INTERMEDIATE  VALUES 
A1,A2,A3,A4,A5,A6/  REAL;  CONSTANTS 
DENOM/  REAL;  DENOMINATOR  OF  THE  EQUATION 


****  PROGRAMMER  /  DAWN  RIGBY  **** 
****  DATE  /  22/06/82  **** 


FUNCTION  RESP(Y) 

IMPLICIT  REAL*8(A-H,  0-Z) 
X  =DABS(Y)  /DSQRT (2 . ODO) 
A1  =  7.05230784D-2 
A2  =  4 . 22820123D-2 


A3  =  9 .27052720-3 
A4  =  1 . 520143D-4 
A5  =  2 . 765672D-4 


A6  =  4 . 30638D-5 

DENOM=  !.♦  A1*X  +  A2*X**2  ♦  A3*X**3  ♦  A4*X**4 


IF (DENOM.LT. 4. 8696750+4) GO  TO  5 
RES  =  0. 

GO  TO  3 

DENOM  =  DENOM **16 
RES  =  . 5/DENOM 
IF(Y-0.)2,2,4 
RESP  =  RES 
GO  TO  6 


RESP  =  1.-  RES 
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A5*X**5  ♦  A6*X**6 
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RETURN 

END 
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****  FUNCTION  PCHISQ  **** 

-  CALLED  TO  EVALUATE  THE  PROBABILITY  FOR  EXCEEDING  CHI  - 

SQUARE,  GIVEN  THE  COMPARISON  VALUE  OF  REDUCED  CHI  SQUARE 
(CHISQR)  AND  THE  NUMBER  OF  DEGREES  OF  FREEDOM  (NFREE) .  THE 
CALCULATION  IS  ONLY  AN  APPROXIMATION  FOR  NFREE  ODD  AND  CHI 
SQUARE  GREATER  THAN  50. 


****  USAGE  **** 

-  RESULT  =  PCHISQ(CHISQR, NFREE)  - 

****  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  *♦** 

-  FUNCTION  GAMMS(X)  - 

CALCULATES  THE  GAMMA  FUNCTION  FOR  INTEGERS  AND  HALF-INTEGERS 

****  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES 

-  PARAMETERS  PASSED  FROM  THE  MAINLINE  - 

CHISQR/  REAL;  COMPARISON  VALUE  OF  REDUCED  CHI  SQUARE 
NFREE/  INTEGER;  NUMBER  OF  DEGREES  OF  FREEDOM 

-  VARIABLES  - 

Z,  TERM,  SUM,  PWR/  REAL;  INTERMEDIATE  VALUES 

NEVEN/  INTEGER;  INTERMEDIATE  VALUE 

FREE  /  REAL;  NUMBER  OF  DEGREES  OF  FREEDOM 

1/  INTEGER;  LOOP  COUNTER 

FI/  REAL;  LOOP  COUNTER 

IMAX/  INTEGER,  MAXIMUM  LOOP  COUNT 


i 

0041 

C 

****  PROGRAMMER  /  DAWN  RIGBY  ** 

0042 

C 

♦  DATE  /  22/06/82  **** 

0043 

C 

0044 

C 

£ 

0045 

FUNCTION  PCHISQ(CHISQR, NFREE) 

0046 

IMPLICIT  REAL*8(A-H,  0-Z) 

% 

0047 

IF  (NFREE) 12 , 12,14 

>7 

0048 

12 

wr i te ( i outdev ,11) 

0049 

11 

FORMAT (3X, 'SUBROUTINE  RESP  --’) 

0050 

wr i te  ( i outdev , 13) 

V, 

0051 

13 

FORMAT (3X, ’**  ERROR  ---  NFREE  ( 

0052 

1  STOP  **’) 

0053 

GO  TO  59 

ra 

0054 

14 

FREE  =  NFREE 

tv 

0055 

Z=  CHISQR*  FREE/2 . 

tv 

0056 

NEVEN=  2  *  (NF  REE/2) 

IF (NFREE  -  NEVLN)?1 ,21 ,41 
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0080 
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56 
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57 
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0087 

59 
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60 

0089 

UNCLASSIFIED 


-  NUMBER  OF  DEGREES  OF  FREEDOM  IS  EVEN  - 

IMAX  =  NFREE/2 
TERM  =  1. 

SUM  =  0. 

DO  34  1=1, IMAX 
FI=I 

SUM=  SUM  ♦  TERM 
TERM  =  TERM  ♦  Z/FI 
PCHISQ  =  SUM*  DEXP(-Z) 

GO  TO  60 

-  NUMBER  OF  DEGREES  OF  FREEDOM  IS  ODD  - 

IF(Z-25.)44,44,42 
Z=  CHISQR*  (FREE-1. )/2. 

GO  TO  21 
PWR  =  FREE/2. 

TERM  =  1. 

SUM  =  TERM/PWR 
DO  56  1=1,1000 
FI=I 

TERM  =  -TERM  *  Z/FI 

SUM  =  SUM  ♦  TERM/ (PWR  ♦  FI) 

IF (DABS (TERM/SUM) - . 00001 ) 57 , 57 , 56 
CONTINUE 

PCHISQ  =  1.0  -  (Z**PWR)  *  SUM/GAMMS (PWR) 

GO  TO  60 
STOP 
RETURN 
END 
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****  FUNCTION  GAMMS  **** 

-  USED  TO  CALCULATE  THE  GAMMA  FUNCTION  FOR  INTEGERS  AND 

HALF-INTEGERS 


****  USAGE  **** 


-  RESULT  =  GAMMS(X)  - 


****  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  **** 

-  FUNCTION  FACTOR (N)  - 

CALCULATES  N  FACTORIAL  FOR  INTEGERS 


*♦**  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES  **** 

-  PARAMETERS  PASSED  FROM  THE  FUNCTION  PCHISO  - 

X/  REAL;  REAL  FORM  OF  INTEGER  OR  HALF-INTEGER  ARGUMENT 
-  VARIABLES  - 


N/  INTEGER;  INTERMEDIATE  VALUE 
XN,SUM/  REAL;  INTERMEDIATE  VALUES 
PROD/  REAL;  CONSTANT 
1/  INTEGER;  LOOP  COUNTER 
FI/  REAL;  LOOP  COUNTER 


****  PROGRAMMER  /  DAWN  RIGBY  **** 
****  DATE  /  22/06/82  ♦*** 


FUNCTION  GAMMS(X) 

IMPLICIT  REAL*8(A-H,  0-Z) 

-  INTEGERIZE  ARGUMENT  - 

N=X- . 25 
XN=N 

IF (X-XN- .75) 31, 31, 21 

-  ARGUMENT  IS  INTEGER  - 

GAMMS  =  FACTOR (N) 

GO  TO  60 

-  ARGUMENT  IS  HALF-INTEGER  - 


PROD  =  1.77245385 
IF(N) 44, 44,33 
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GAMMS 

0058 

33 

IF(N-10)41,41,51 

0059 

41 

DO  43  1=1, N 

0060 

FI  =  I 

0061 

43 

PROD  =  PROD  *  (FI- .5) 

0062 

44 

GAMMS  =  PROD 

0063 

GO  TO  60 

0064 

51 

SUM  =  0. 

0065 

DO  54  1=11, N 

0066 

FI  =  I 

0067 

54 

SUM  =  SUM  +DL0G(FI-.5) 

0068 

GAMMS  =  PROD  *  639383.8623  *DEXP(SUM) 

0069 

60 

RETURN 

0070 

END 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 


****  FUNCTION  FACTOR  **** 

-  CALLED  TO  CALCULATE  THE  FACTORIAL  FUNCTION  FOR  INTEGERS  - 


****  USAGE  **** 

-  RESULT  =  FACTOR (N) 


>**  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED 


-  NONE  — 


****  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES  **** 

-  PARAMTERS  PASSED  FROM  THE  FUNCTION  GAMMS  - 

N/  INTEGER;  INTEGER  ARGUMENT  FOR  CALCULATIONS 

-  VARIABLES  - 

1/  INTEGER;  LOOP  COUNTER 
FI/  REAL;  LOOP  COUNTER 
SUM/  REAL;  INTERMEDIATE  VALUE 


****  PROGRAMMER  /  DAWN  RIGBY  **** 
****  DATE  /  22/06/82  **** 


FUNCTION  FACTOR (N) 

IMPLICIT  REAL*8(A-H,  0-Z) 
FACTOR  =  1. 

IF(N-1)40,40,13 
IF(N-10)21 ,21,31 
DO  23  1=2, N 
FI  =  I 

FACTOR  =  FACTOR  *  FI 
GO  TO  40 
SUM  =  0. 

DO  34  1=11, N 
FI  =  I 

SUM  =  SUM  ♦  DLOG(FI) 

FACTOR  =  3628800. *DEXP (SUM) 

RETURN 

END 
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0001  C 

0002  C 

0003  C  **♦*  FUNCTION  TPROB  ♦*** 

0004  C 

0005  C  -  THIS  FUNCTION  FINDS  THE  VALUE  OF  »T*  FROM  A  *T»  - 

0006  C  DISTRIBUTION  OF  'NFREE'  DEGREES  OF  FREEDOM,  WHICH  IS  EXCEEDED 

0007  C  WITH  PROBABILITY  ALPHA  (2  SIDED  TEST) .  A  TABLE  IS  USED  AS  A 

0008  C  DATA  ARRAY  AND  INTERPOLATION  IS  OFTEN  NECESSARY  TO  GET  THE 

0009  C  CORRECT  VALUE  OF  TPROB.  THE  TABULATED  DATA  WAS  COMPUTED  BY 

0010  C  MAXINE  MERRINGTON,  FROM  "TABLES  OF  PERCENTAGE  POINTS  OF  THE 

0011  C  INCOMPLETE  BETA  FUNCTI0N\BI0METRIKA,VGL.32  (1941), PAGES  168- 

0012  C  181;  BY  CATHERINE  M.  THOMPSON. 

0013  C 

0014  C 

0015  C  ****  USAGE  **** 

0016  C 

0017  C  -  RESULT  =  TPROB (FREE .ALPHA)  - 

0018  C 

0019  C 

0020  C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  CALLED  **** 

0021  C 

0022  C  - NONE - 

0023  C 

0024  C 

0025  C  ♦***  DESCRIPTION  OF  PARAMETERS  AND  VARIABLES  *♦** 

0026  C 

0027  C  -  PARAMETERS  PASSED  FROM  THE  MAINLINE  - 

0028  C 

0029  C  FREE/  REAL;  NUMBER  OF  DEGREES  OF  FREEDOM 

0030  C  ALPHA/  REAL;  CONFIDENCE  LEVEL  OF  THE  FIDUCIAL  LIMITS 

0031  C 

0032  C  - VARIABLES - 

0033  C 

0034  C  T /  REAL;  ARRAY  FOR  T  DISTRIBUTION  TABLE 

0035  C  A/  REAL;  1-ALPHA,  FOR  READING  OF  THE  TABLE  AND  CALCULATIONS 

0036  C  I,J /  INTEGER;  LOOP  COUNTERS  AND  ARRAY  POINTERS 

0037  C  AK,AJK,H,YK,YJK,DY,TI,TIA,FR/  REAL;  INTERMEDIATE  VALUES 

0038  C  FARRAY/  REAL;  ARRAY  FOR  INTERPOLATED  FREE  VALUES  WHEN  BOTH 

0039  C  FREE  AND  A  MUST  BE  INTERPOLATED 

0040  C 

0041  C 

0042  C  ****  PROGRAMMER  /  DAWN  RIGBY  **♦* 

0043  C  ****  DATE  /  22/06/82  **** 

0044  C  ****  Updated  /  John  McFee  **** 

0045  C  ****  DATE  /  26/02/86  **** 

0046  C 

0047  C 

0048  FUNCTION  TPROB (FREE , ALPHA) 

0049  IMPLICIT  REAL*8(A-H,  0-Z) 

0050  DIMENSION  T (35, 8), FARRAY (7) 

0051  data  itabledev/4/ 

0052  READ  ( i  tab  I  edev ,  1)  (  (T  ( J ,  I) ,  1=1 , 8) ,  J=1 , 35) 

0053  1  FORMAT (8F10. 5) 

0054  A  =  l.ODO  -  ALPHA 

0055  FR=FREE 

0056  C 

0057  C  -  CHECK  FOR  OUT  OF  RANGE  ’A’  OR  ’FREE’  VALUES  - 

UNCLASSIFIED 
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TPROB 

0058 

C 

0059 

0060 

0061 

0062 

3 

0063 

2 

0064 

0065 

9 

0066 

0067 

0068 

0069 

5 

0070 

0071 

11 

0072 

0073 

7 

0074 

0075 

0076 

13 

0077 

0078 

19 

0079 

0080 

0081 

0082 

15 

0083 

0084 

21 

0085 

0086 

17 

0087 

0088 

23 

0089 

0090 

C 

0091 

C 

0092 

C 

0093 

0094 

0095 

27 

0096 

C 

0097 

C 

0098 

C 

0099 

28 

0100 

0101 

31 

0102 

35 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

C 

0112 

C 

0113 

C 

0114 

25 
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IF (A. LT. 0.005) GO  TO  3 
IF (A. GT. 0.50) GO  TO  5 
GO  TO  7 

write(ioutdev,2) 

FORMAT (//3X, ’SUBROUTINE  TPROB — ’) 
wri te(ioutdev,9) 

FORMAT (3X, ’ALPHA  IS  GREATER  THAN  0.995,  THEREFORE,  ASSUME  THAT  ALP 
1HA  IS  EQUAL  TO  0.995’) 

A=0 . 005 
GO  TO  7 

write(ioutdev,2) 

write(ioutdev,ll) 

FORMAT(’l’, ’ALPHA  IS  LESS  THAN  0.50  —  ERROR  ’) 

GO  TO  34 

IF (FR . GT . 120 . 0) GO  TO  13 
IF(FR.LT.1.0)GO  TO  15 
GO  TO  17 
write('ioutdev  ,2) 
wr ite(ioutdev,19) 

FORMAT (3X, ’THE  NUMBER  OF  DEGREES  OF  FREEDOM  IS  GREATER  THAN  120.0; 
l’/10X, ’THEREFORE  ASSUMED  TO  BE  INFINITE’) 

FR  =  900. 

GO  TO  17 
write(ioutdev,2) 
wr ite(ioutdev ,21) 

FORMAT(’l ’ , 'FREE  IS  LESS  THAN  1.0  —  ERROR  ’) 

GO  TO  34 

IF(FR.LE.30.0)GO  TO  23 
GO  TO  25 
DO  27  1=2,8 

IF(T (1 , I) ,NE.A)GO  TO  27 

-  FREE  AND  A  ARE  BOTH  ON  TABLE  - 

TPROB  =  T (FR+1 , I) 

GO  TO  33 
CONTINUE 

-  A  IS  NOT  ON  TABLE  - 

DO  31  1=2,8 

IF  (T  (1,1) . LT . A) GO  TO  35 

CONTINUE 

AK=T (1 , I) 

1=1-1 

AJK=T (1,1) 

H=AJK-AK 

YJK=T(FR*1,I) 

YK=T(FR*1,I*1) 

DY=YJK-YK 

TPROB  =  (DY  *  (A-AK)/H)  +YK 
GO  TO  33 

-  FREE  IS  GREATER  THAN  30.0  - 

DO  37  J=31 ,35 
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0115 

IF(T(J,1).NE.FR)C0  TO  37 

0116 

DO  75  1=2,8 

0117 

TI=T(1,I) 

0118 

TIA  =  TI  -  A 

0119 

IF (TIA . LT . 1 . OD-13)  GO  TO  76 

0120 

75 

CONTINUE 

0121 

GO  TO  28 

0122 

C 

0123 

C 

-  FREE  >  30.0;  FREE  AND  A  B 

0124 

C 

0125 

76 

TPROB  =  T(J,I) 

0126 

GO  TO  33 

0127 

37 

CONTINUE 

0128 

39 

DO  41  1=2,8 

0129 

IF (T (1,1) .NE.AJGO  TO  41 

0130 

C 

0131 

C 

-  FREE  IS  NOT  ON  TABLE  - 

0132 

C 

0133 

DO  42  J=31,35 

IF(T(J,1) .GT.FR)GO  TO  45 

0134 

0135 

42 

CONTINUE 

0136 

41 

CONTINUE 

0137 

GO  TO  43 

0138 

45 

AK=T ( J , 1) 

0139 

J=J-1 

0140 

AJK=T ( J , 1) 

0141 

H=AJK-AK 

0142 

YJK=T ( J , I) 

0143 

YK=T (J*1,I) 

0144 

DY=YJK-YK 

0145 

TPROB  =  (DY  ♦  (FR-AK) /H)  +YK 

0146 

GO  TO  33 

0147 

C 

0148 

C 

-  NEITHER  FREE  NOR  A  ARE  ON 

0149 

C 

0150 

43 

DO  47  J=31,35 

0151 

IF(T (J, 1) .GT.FR)GO  TO  49 

0152 

47 

CONTINUE 

0153 

49 

AK=T (J, 1) 

0154 

J=J-1 

0155 

AJK=T(J,1) 

0156 

H=AJK-AK 

0157 

DO  51  1=2,8 

0158 

YJK=T ( J , I) 

0159 

YK=T ( J>1 , I) 

0160 

DY=YJK-YK 

0161 

TPROB  =  (DY  *  (FR-AK) /H)  ♦  YK 

0162 

FARRAY (1-1) =TPROB 

0163 

51 

CONTINUE 

0164 

DO  53  1=2,8 

0165 

IF (T (1,1) .LT.A) GO  TO  55 

0166 

53 

CONTINUE 

0167 

55 

AK=T(1,I) 

0168 

1=1-1 

0169 

AJK=T(1,I) 

0170 

H=AJK-AK 

0171 

YK  =  FARRAY (I- 1) 
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0172 

0173 

0174 

0175 

0176  34 
0177  33 


YJK  =  F ARRAY (1-2) 

DY=YJK-YK 

TPROB  =  (DY  *  (A-AK) /H)  ♦  YK 

GO  TO  33 

STOP 

RETURN 


A5/23 

26-Feb-1986  16:03:43  VAX  FO 
26-Feb-1986  14:18:08  DISKSM 


UNCLASSIFIED 
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I 


OPERATIONAL  DESCRIPTION  OF 
RELATIVE  POTENCY  PROGRAM  RPOTC 
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RELATIVE  POTEMCv  PROGRAM  (P.POTC) 

Revised  Edition  —  Account  BI000200, STATS, OPES  — 

Hawn  Rigby  —  June,  1983 


INPUT: 

The  RPOTC  program  has  now  been  transferred  to  the  Honeywell  CP-6 
computer  system,  and  some  changes  in  input  are  therefore  required.  The  data 
is  typed  interactively  (excluding  the  chisq  and  t  table)  and  read  directly  by 
the  command  file.  Each  value  is  asked  for  separately,  so  that  carriage 
returns  separate  the  data.  The  values  are  as  follows. 

ISLCT  -  select  value.  If  ISLCT=1,  6  sets  of  numbers  must  follow  (limit, 
lowest  datum,  and  highest  datum  for  each  set).  This  indicates  that  limits 
other  than  %%  (ISLCT=Q)  are  to  be  used.  If  a  -1  is  typed,  execution 
terminates. 

LIMIT  -  percentage  limits  to  be  used  (integer) 

I  Min  -  number  of  comparisons  to  be  made  in  the  set 
I  -  lowest  datum 
J  -  highest  datum 
I CODE  -  identification  code 

CHISQ ,SNV,XBAR ,SXX ,SXY .HFREE ,YBAR ,SYY  -  results  of  the  PROBIT  analysis 
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EXECUTION 

At  this  time,  the  mainline  RD0TC  program  is  in  the  file  KEEP_RPOTC,  and 
the  command  file  to  run  this  is  called  RPOTC.  These  files  MUST  NOT  be 
deleted  under  any  circumstances .  The  compilation  and  linage  files  are 

called  RPOTC _ OU  and  R°0TC_RU.  These  should  also  be  left  alone,  but  can  be 

replaced  by  typing  the  following: 

FORTRAN  KEEP_RPOrC  OVER  RP0TC_0ll 

LINK  RPOTC  OU  OVER  RPOTC  RU 


To  execute  this  program,  type: 

XEO  RPOTC 

The  program  and  additional  files  are  all  in  the  acount 
3 1000200, STATS, ORES. 


OUTPUT 


The  output  for  this  program  has  not  been  changed,  and  is  sent  to  both 
the  terminal  and  the  line  printer  immediately  after  execution. 
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TABLE  1 

Initial  Data  for  Program  PROBIT 


ID  -  array  for  the  identification  label  (to  begin  each  page) 
format:  alphabetic;  maximum  of  50  characters 

NT  -  number  of  trials,  maximum  of  100  (n) 

MAXCNT  -  maximum  loop  count  for  iterative  regression  (default  »  50) 
format:  integer;  NT,  MAXCNT  are  entered  on  the  same  line. 

IYNSW  -  switch  for  output  of  intermediate  results 
YES  -  IY,  output  requested 
NO  -  IN;  no  output  requested 
format:  alphanumeric;  maximum  of  50  characters 

TN  -  number  of  animals  in  the  trial  (n^) 

XIN  -  dose  given  in  the  trial  ( x ! ) 

R  -  number  of  responses  in  the  trial,  maximum  of  100  (r..) 

format:  real;  TN,  XIN,  R  are  entered  on  the  same  line,  1 
trial /line 

ALPHA  -  confidence  level  of  fiducial  limits  (default  =  0.95),  { a/100) 
format:  real 

C  -  fraction  of  natural  mortality  (default  *  0.0)  (PQ) 

DELTA  -  convergence  criterion  (default  =  0.001)  (s) 

format:  real;  both  C  and  DELTA  typed  on  the  same  line 
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TABLE  2 


Quantities  Output  by  Program  PROBIT 


I 

J 

X 

RP 

RS 

YK 

YR 

SIGY 

W 

B 

SIGMAB 

FLY 

FUL 

CHISQ 

NFREE 

H 


-  real 

-  Real 

-  real 

-  real 

-  real 

-  real 

-  real 

-  real 

-  real 

-  real 


integer;  trial  number  (i) 
integer;  loop  number  (k) 
real . 

real;  empirical  percentage  responses  (P.*10Q) 

1  A 

mean  of  estimated  number  responding  (R^) 
maximum  likelihood  expected  NEDS  (y(x-)) 
probits  (NEDS  +  5) 
standard  error 


weighting  coefficient  (w- 


(k) 


) 

(k) 


) 

(k) 


T 

-  real 

G 

-  real 

PLB 

-  real 

EOLMB 

-  real 

XLMBDA 

-  real 

FUED 

-  real 

FLED 

-  real 

FUX 

-  real 

FLX 

-  real 

APA 

-  real 

current  regression  slope  (a, 
error  in  regression  slope  (S^'*') 
lower  fiducial  limits  of  estimated  number  responding 
upper  fiducial  limits  of  estimated  number  responding 
value  of  chi  square  ( (  v) ) 

-  integer;  number  of  degrees  of  freedom  (  vj) 

-  real;  heterogeneity  factor  (h) 

H  =  1.0  ---  homogeneous 
H  *  1.0  —  deterogeneous 

value  of  T  distribution  (t  ) 

G  factor  (g) 


a 


given  ED  levels  (dose  level  *  100)  (P^l 
tolerance  level  in  dose  units  (ED^) 
tolerance  level  in  dosage  units  ( x x) 
upper  fiducial  limit  for  EDLMB  (f^) 
lower  fiducial  limit  for  EDLMB  (fj  x) 
upper  fiducial  limits  for  XLMBDA  (f^) 
lower  fiducial  limits  for  XLMBDA  (f^) 
percent  confidence  level  (a) 
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TABLE  3 

Data  for  Comparison  Tests  of  Programs  PROBIT  and  PROBT 


Dose 

( Relative  Units) 
Series  1 


Series  2 


Series  3 


Series  4 


Series  5 


Organisms  Responses 


Dose 

(Relative  Units) 
Series  6 


Series  9 


Organisms  Responses 


150. 

20 

19 

140. 

30 

26 

130. 

40 

29 

120. 

50 

27 

no. 

40 

14 

100. 

30 

6 

80. 

20 

0 

Series  7 

499. 

10 

4 

474. 

10 

5 

450. 

20 

6 

427. 

10 

1 

406. 

10 

0 

405. 

10 

1 

364. 

10 

1 

327. 

10 

0 

Series  8 

499. 

10 

6 

474. 

10 

6 

450. 

20 

6 

427. 

10 

1 

406. 

10 

0 

405. 

1U 
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TABLE  5 

Data  for  Comparison  Tests  of  Programs  PROBIT,  S103  and  SAS 


Dose 

(RelaTTve"  Units) 


Organisms  Responses 


Dose 

(Relative  Units) 


Organisms 


Responses 


1 

i 


Series  1 

6374. 

10 

10 

5951. 

10 

8 

5371. 

10 

1 

4975. 

10 

2 

4522. 

10 

2 

4048. 

10 

0 

3674. 

10 

1 

3312. 

10 

0 

Series  2 

5951. 

10 

10 

4975. 

10 

7 

4522. 

10 

7 

4048. 

10 

0 

Series  3 

6374. 

10 

4 

5951. 

10 

1 

5371. 

10 

1 

4975. 

10 

1 

4522. 

10 

0 

4048. 

10 

0 

3674. 

10 

0 

3312. 

10 

0 

Series  4 

6374. 

10 

10 

5951. 

10 

10 

5371. 

10 

3 

4975. 

10 

6 

4522. 

10 

2 

4048. 

10 

0 

3674. 

10 

1 

3312. 

10 

0 

Series  5 


6374. 

10 

10 

5971. 

10 

10 

5371. 

10 

2 

4975. 

10 

3 

4522. 

10 

2 

4048. 

10 

0 

3674. 

10 

1 

3312. 

10 

0 

Series  6 

50. 

10 

9 

40. 

10 

8 

32. 

10 

9 

25. 

10 

0 

Series  7 

650. 

10 

10 

585. 

10 

10 

526.5 

10 

6 

473.9 

10 

2 

426.5 

10 

3 

Series  8 

222. 

10 

10 

200. 

10 

3 

180. 

10 

2 

162. 

10 

0 

146. 

10 

0 

131. 

10 

0 

118. 

10 

0 

106. 

10 

0 

Series  9 

585. 

10 

10 

526.5 

10 

6 

473.9 

10 

2 

426.5 

10 

3 
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